Two spin liquid phases in the spatially anisotropic triangular Heisenberg model 
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The quantum spin- 1/2 antiferromagnetic Heisenberg model on a two dimensional triangular lat- 
tice geometry with spatial anisotropy is relevant to describe materials like CS2CUCI4 and organic 
compounds like k-(ET)2Cu2(CN)3. The strength of the spatial anisotropy can increase quantum 
fluctuations and can destabilize the magnetically ordered state leading to non conventional spin 
liquid phases. In order to understand these intriguing phenomena, quantum Monte Carlo methods 
are used to study this model system as a function of the anisotropic strength, represented by the 
ratio J' /J between the intra-chain nearest neighbor coupling J and the inter-chain one J'. We have 
found evidence of two spin liquid regions. The first one is stable for small values of the coupling 
J' I J < 0.65, and appears gapless and fractionalized, whereas the second one is a more conventional 
spin liquid with a small spin gap and is energetically favored in the region 0.65 < J'/J < 0.8. We 
have also shown that in both spin liquid phases there is no evidence of broken translation symmetry 
with dimer or spin-Peirls order or any broken spatial reflection symmetry of the lattice. The various 
phases are in good agreement with the experimental findings, thus supporting the existence of spin 
liquid phases in two dimensional quantum spin- 1/2 systems. 

PACS numbers: 71.10.-w,71.10.Pm,75.10.-b,75.40.Mg 



I. INTRODUCTION 

Since the pioneering work by Anderson and Fazekas 1 
the spin- 1/2 antiferromagnetic Heisenberg model on the 
triangular lattice has been considered one of the most 
promising candidates for a spin liquid phase in a frus- 
trated antiferromagnet. However several numerical stud- 
ies 2 ~ 5 have all consistently confirmed that in the isotropic 
triangular lattice the classical magnetically ordered state 
appears stable. Nevertheless the ordered moment is 
found considerably smaller than the classical value, 4,5 
suggesting that the model is very close to a quantum 
critical point, 6 namely to a phase where the long range 
antiferromagnetic order is completely destroyed. This 
picture is supported by the recently established result 
that in the quantum dimer model on the triangular lat- 
tice geometry 7 a spin liquid phase is stable, a result that 
is particularly important because for instance the same 
model on the square lattice displays only non spin liquid 
phases with broken translation symmetry. 8-10 

Recent experiments on two different materials 11 ' 12 
have renewed the interest in the spin- 1/2 antiferromag- 
netic Heisenberg model on the triangular lattice de- 
scribed by the following Hamiltonian: 13 



chains, implying spin fractionalization and no antiferro- 
magnetic order. In this limit the spin one excitations 
should form a broad two-spinon continuum of states as 
predicted theoretically, 14 and indeed several experiments 
by inelastic neutron scattering have revealed this non 
trivial spin dynamics. 15 
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where S\ is a spin 1/2 located at site i on the triangu- 
lar lattice, indicates nearest neighbor sites 
along the chain (between different chains), and the cor- 
responding antiferromagnetic couplings are denoted by J 
and J' (see Fig. 1). 

Clearly the anisotropy increases the quantum fluctu- 
ations in this model as for J' = the Hamiltonian H 
reduces to a system of uncoupled one dimensional (ID) 



FIG. 1: Antiferromagnetic Heisenberg models on the trian- 
gular lattice studied. A spin-1/2 is located at each dot. 

?i = (1,0) and T2 = (5, 2 ) denote the primitive transna- 
tional vectors. A lattice constant is set to be one. 



In a series of experiments by Coldea, et al., 11 the low 
energy spin dynamics of CS2CUCI4 have been studied sys- 
tematically and found to be described essentially by the 
model Hamiltonian (1) with anisotropic exchange inter- 
action J I J' ~ 3.0. In their experiments, 11 there are two 
facts supporting the existence of a spin liquid phase in 
this system: i) a spin spiral phase appears at tempera- 
tures below Tjv — 0.6 K, which is about the same order 
of magnitude as the inter-plane coupling J" / J ~ 0.045 in 
the third spatial direction, and therefore sizably smaller 
than the two dimensional (2D) couplings J (~ 0.37 meV) 
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and J'. In fact, true long range order in 2D is possible 
only at zero temperature and a finite Tm is due to the in- 
ter plane coupling J" which allows to cut-off the logarith- 
mically divergent quantum fluctuations in 2D. A finite 
Tjv is therefore expected to be of order ~ J'/ log(J'/ J"), 
which appears slightly off from the observed Tm — 0.6 K. 
ii) At temperatures larger than Tm (or by applying an 
external magnetic field), the inelastic neutron scattering 
experiments 11 found that the line shape of the spectrum 
consists of a broad continuum, which is in contrast to 
the expected behavior of a magnetically ordered state, 
but is instead similar to the broad two-spinon continuum 
expected in the ID systems mentioned above. Therefore 
these experiments suggest that spin fractionalization can 
be realized also in this 2D system for temperatures higher 
than the corresponding 3D transition temperature Tm- 

Another spectacular experiments on spin-1/2 quasi 
2D triangular systems have been recently reported 
by Shimizu et al. 12 Here two organic materials k- 
(ET) 2 Cu 2 (CN) 3 and k-(ET) 2 Cu 2 [N(CN) 2 ]C1 (where ET 
denotes bis(ethylenedithio)-tetrathiafulvalencc) are syn- 
thesized, corresponding essentially to two different values 
of J' /J ~ 0.8-0.9 and J' /J ~ 1.8, respectively 12 While 
the latter material shows commensurate spin ordering at 
Tm = 27 K, for the former material no magnetic order 
is observed even down to the milli-Kclvin region (~ 32 
mK) regardless the fact that the estimated value of J is 
about 250K. This observation strongly indicates possible 
realization of a spin liquid state for k-(ET) 2 Cu 2 (CN)3. 12 

Motivated by these experiments, we shall consider here 
the spin-1/2 antiferromagnetic Heisenberg model on the 
triangular lattice with spatial anisotropy described by 
Eq. (I). In particular, we report a detailed and system- 
atic quantum Monte Carlo (QMC) study of the ground 
state as well as the low-lying excitations of this model 
system as a function of J' /J. Using both quantum vari- 
ational Monte Carlo (VMC) method and Green function 
Monte Carlo method with an improved extension of the 
fixed node approximation (named effective Hamiltonian 
approach in this paper), it is found that there exist two 
spin liquid regions in the phase diagram of this model 
by varying the ratio J'/ J for J' / J < 0.8. The first one 
is stable for small values of the coupling J' / J < 0.65, 
and appears to show gapless, fractionalized fermionic ex- 
citations. The second one is energetically favored in the 
region 0.65 < J'/ J < 0.8, and is a more conventional spin 
liquid with a small spin gap in the excitation spectrum, 
the same type of spin liquid phase realized in the quan- 
tum dimer model in the isotropic triangular geometry. 7 It 
is argued that the two experimental observations of spin 
liquid like behaviors for Cs 2 CuCl 4 and k-(ET) 2 Cu 2 (CN) 3 
mentioned above should correspond to these two different 
spin liquid phases of the model, respectively. 

The paper is organized as follows. In Sec. II, we 
first introduce the variational wave functions consid- 
ered (Sec. II A), and describe our original optimization 
method to obtain the minimum energy variational wave 
function containing a large number of variational param- 



eters (Sec. II B). In order to systematically correct this 
variational ansatz, an effective Hamiltonian approach is 
introduced in Sec. II C along the line of the well estab- 
lished diffusion Monte Carlo technique, 16 allowing, for 
continuous systems, to achieve the best variational wave 
function with the same phases of the chosen variational 
state. As explained in Sec. II D, within this approach, 
it is also possible to calculate the low-lying excitation 
spectrum, which can be compared directly with dynam- 
ical experimental measurements. In Sec. Ill, all our nu- 
merical results are reported for the spin-1/2 antiferro- 
magnetic Heisenberg model on the triangular lattice as 
a function of J' / J , including the decoupled chains case 
(J 7 = 0). Finally our conclusions and remarks are pre- 
sented in Sec. IV. The paper is also supplemented by 
several important appendices for the detailed explana- 
tion of the methods and the wave functions used, which 
should be very useful for reproducing our results or ex- 
tending our approach to other model systems. A part of 
this work have been reported briefly as a short commu- 
nication. 17 



II. NUMERICAL METHOD 

A. Variational wave functions: projected BCS 
states 

The variational wave function considered in this study 
is the so-called projected BCS state defined by 



|p-BCS) =7> G exp 



E 

i<j 



fij ( C !,T C L + c J,T c a) 



|o>, 



(2) 

where cJ CT (cj ;Cr ) is an electron creation (annihilation) op- 
erator at site i with spin cr(=|,|), |0) is the vacuum 
state with no electrons, the function fij is the so-called 
pairing function that contains all variational freedom of 
the |p-BCS), and is determined by the minimum energy 
condition, whereas Vq is the usual Gutzwillcr projection 
operator onto the subspace of singly occupied sites, im- 
plying that the total number of electrons N is equal to 
the number of sites L. The pairing function fij of this 
projected BCS state can be parameterized using a BCS 
Hamiltonian: 



(3) 

Here tij and Aij as well as the chemical potential 
U,i = — A* (assumed uniform) can be considered varia- 
tional parameters, which implicitly determine the pair- 
ing function fij corresponding to the ground state (GS) 
of Hbcs- Here i,j (= 1,2, ■■ ■ , L) label the sites of the 
lattice (see Fig. 1), i.e., Fj = i\T X + i 2 r 2 , in the lexico- 
graphic order, so that the condition i < j in Eq. (2) is 
meaningful in any spatial dimension. 
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When tij and Ajj depend only on I = fi — fj, i.e., 

tij = tf and Aij — Af, respectively, Hbcs can be de- 
scribed more compactly by 

#bcs = ^(e k -M)4,a c k^+I](Ak4 T cL ka +h.c). (4) 



k.cr 



Here 
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(5) 



(6) 



(7) 



For a singlet pairing, A it j — A^j, and thus A k = A_ k . 
In this case the projected BCS wave function defined by 
Eq. (2) reads 



|p-BCS) =7> g |BCS), 



(8) 



opposed to the ID case, may differ from zero in the trian- 
gular lattice case. As will be discussed in Sec. IIIC, we 
found that the variational parameters t? 2 and t? 3 are ir- 
relevant and the energetically favorable symmetry of A k 
is A x . 17 

As shown in Sec. Ill, the projected BCS state de- 
scribed above is a very good variational state for J'/ J < 
0.7. However, close to the isotropic limit J' /J ~ 1.0, 
the translation invariant ansatz state, also previously at- 
tempted in the presence of hole doping, 19 ' 20 is not very 
accurate when compared with the exact diagonalization 
results possible on the 6x6 cluster. 21 In this region of 
J'/ J, as shown in Sec. HID, we have found that it is 
more convenient to consider a BCS Hamiltonian defined 
on a (2 x 1) unit cell [cf. Eq. (55)] for a much better 
variational wave function. As shown in App. C, by using 
the projected BCS state thereby constructed, it is then 
possible to represent the well known short range reso- 
nant valence bond (RVB) wave function, with a particu- 
lar choice of the variational parameters. The RVB state 
is a good variational ansatz for the isotropic triangular 
lattice 3 and represents a very convenient initial guess for 
defining, within the present |p-BCS) framework, an ac- 
curate variational state in the nearly isotropic triangular 
lattice. 



with |BCS) being the ground state of H BCS given by 
Eq. (4): 



Minimization method 



| BCS) = cxp 



EM: 



T C k,i 
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where / k = v k /u k = A k / (£ k + E k ), u k 



1 + «!S- 



and 



(10) 



with £k = £k — A*- 

At the variational level, both £ k and A k have to be 
parameterized in order to minimize the variational en- 
ergy of the |p-BCS) wave function. The most rele- 
vant parameters for lowering the energy are the short 
range terms, and we have chosen to expand = 

J2t=i 2*t„ cos (k ■ T^j , where tf v are variational parame- 
ters and t v are the nearest neighbor vectors (T3 = T2— t±). 
Analogously, the gap function A ; - is truncated up to the 
third nearest distance along the chain (ri) direction. This 
is also because, as will be discussed in Sec. IIIB, for the 
ID spin-1/2 antiferromagnetic Heisenberg model Hid, 
inclusion of the parameter A 3 ^ is known to be crucial 
for this type of projected BCS states to represent almost 
exactly the ground state of Hid. 18 For the present 2D 
system described by Eq. (1), which preserve Ci v symme- 
try for J' 7^ J, the projected BCS state |p-BCS) thus 
contains ten independent variational parameters in A k 
(see, e.g., Tab. I) and the chemical potential \i which, as 



In order to evaluate the optimal variational parameters 
that minimize the energy expectation value, 



(mm 



(ii) 



we follow the method, recently introduced for calcula- 
tions of electronic structure, 22 which will be described in 
some detail in the following. According to this method, in 
order to reach the minimum energy E(^f) in a stable and 
efficient way, the logarithmic derivative O k (x) of the wave 
function tyr ak y(x) — {x\^i ak y} with respect to each vari- 
ational parameter a k = t^j and/or Aij (k = 1, 2, . . . ,p) 
has to be evaluated on a given iV-electron configuration 
x = {ri, . . . , r^} of the projected Hilbert space with one 
electron per site, namely, 



d 

Ok(x) = —\^ {cth} {x). 



(12) 



In fact, for infinitesimal changes of these variational pa- 



rameters, Ol k 



= a k + Seek, the corresponding change 



of the wave function reads 



l + j20 k (x)-6a k + 0(5a 2 k ) 



k=l 



(13) 

In order to simplify the notations, here we introduce the 
operator O k corresponding to O k (x) which is defined by 



(x\O k \x'} = O k (x) ■ 5 X 



(14) 
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Using this operator form, the above equation is more 
compactly written as 



Sa k O k 



fc=i 



l*W}> 



(15) 



valid up to 0(<5a|). 

For the minimization method described below, it is 
important to evaluate numerically the value O k (x) cor- 
responding to a given real space configuration of elec- 
trons x satisfying the constraint of no doubly occupied 
sites. To this purpose, we have to recall that the varia- 
tional parameters Ajj and tij are explicitly defined in 
Hbcs [Eq. (3)], but only implicitly in the wave function 
itself, defined as the ground state of this unprojected BCS 
Hamiltonian |BCS) [Eq. (2)]. Thus, in order to evaluate 
the logarithmic derivatives of the wave function l^^) 



with respect to a variational parameter a k , we apply sim- 
ple perturbation theory to -£/bcs an d calculate the per- 
turbed state \^{a k +Sa k }), within linear order in Sa k . It 
is then possible to compute (x\^ i ak+ g ak y) using simple 
algebra and within the same accuracy O(Sal). After 
recasting the calculation by using appropriate matrix- 
matrix operations, the evaluation of O k (x) is possible in 
an efficient way, using only ~ (2L) 2 operations for each 
variational parameter a k . More details of the method 
are found in App. A. 



The minimization method used here is similar to the 
standard and well known steepest descent method, where 
the expectation value of the energy ^(*{ Qfc }) is opti- 
mized by iteratively changing the parameters a k (k = 
1, • • • ,p) according to the corresponding derivatives of 
the energy (generalized forces): 



fk 





(*{a»}l 


'6 k H + HO k 


l*w>) 


da k 


<*{a»}l*{a»}) 





+ 2E(V {ak} ) 



(16) 



namely, 



a k 



a' k = a k + f k ■ St. 



(17) 



where St in the standard steepest descent method is de- 
termined at each iteration by minimizing the energy ex- 
pectation value. 23 The above expressions for the forces f k 
can be computed statistically by appropriately averaging 
the local energy 

er(x)- <*w|g|*> 

and O k (x) over a set of configurations {x;}, I = 
1,2, ■■■ , M distributed according to the square of the 

I 1 2 

wave function, (a;|^{ ajk }) , generated with a standard 
variational Monte Carlo (VMC) scheme, namely, 



M 



fk 

O k 



~M °k( x i) e L{xi) + 20 k e L 
i=i 

1 M 



M 



1 M 

i=i 



(18) 
(19) 
(20) 



and thus f k can be computed efficiently once O k (x) and 
e^(x) are evaluated for any sampled configuration. Here 
we assumed that all quantities are real. However an ex- 
tension to the complex case is straightforward. 



Now, for simplicity, we assume that St is positive and 
small enough. Indeed the variation of the total energy 



AE = E{9 Wk} )-E(9 {ak} ) 
= -J2f k -Sa k +0(Sa 2 k ) 



(21) 



k=l 



at each step is easily shown to be negative for small 
enough St because in this limit 



AE=-Stj2f k +0(St 2 ). 



(22) 



k=l 



Thus the steepest descent method certainly converges to 
the minimum of the energy when all the forces vanish. 

Let us now generalize the steepest descent method by 
slightly modifying the basic iteration (17) with a suitably 
chosen positive definite matrix s: 



a k -> a' k = a k + St^s k jfi. 



(23) 



Again, using the analogy with the steepest descent 
method, convergence to the energy minimum is reached 
when the value of St is sufficiently small and is kept pos- 
itive constant for each iteration. In fact, similarly to 
the steepest descent method, the energy variation corre- 
sponding to a small change of the parameters is: 



AE = -Stj2j2s^fkfi + 0(St 2 ). 



(24) 



fe=i 



■5 



and is always negative for small enough St, unless the 
minimum condition of f k = is reached and the varia- 
tional parameters no longer change because Sa k = 0. It 
should be noted here that the steepest descent method is 



J 



a special case with the matrix s = 1 (unit matrix). 

A more convenient choice for the matrix Sj k is given 
by 22 



Sj.k — 



(*{a fc }|OjOfc|*{a»}> (*{a fc }|Oj|*{a»}> (*{a»}|6fc|*{a fc }> 



<*{a*}l*{a*}> 



(25) 



r 



This matrix can be efficiently evaluated statistically, and 
similarly to the forces fk , can be obtained within a VMC 
scheme once Ok(x) are computed for a set of configura- 
tions {xi}, I — 1, 2, ■ ■ ■ ,M distributed according to the 

I 1 2 

wave function squared (x|4 r { Qfc }) , namely, 



1 M 



(26) 



i=i 



where Ok and Oj are defined in Eq. (19). For what- 
soever choice of the M configurations {xi}, this matrix 
remain positive definite regardless of the statistical noise 
(at most has vanishing eigenvalues), because it is explic- 
itly written in Eq. (26) as an overlap matrix in 1Z M be- 
tween the p vectors Ok(xi) — Ok (k = 1, • • • ,p). In this 
paper, this generalized method defined by Eq. (23) with 
the matrix s given by Eq. (25) will be called Stochastic 
Reconfiguration (SR) optimization method. It should be 
noted here that other convenient types of positive def- 
inite matrix have been recently proposed. 24,25 However 
for the present purpose the improvement does not ap- 
pear important. 

Let us next discuss why the choice Eq. (25) for the ma- 
trix s in the SR scheme [Eq. (23)] is particularly simple 
and convenient compared to the simplest steepest descent 
method. For a stable iterative method for the energy op- 
timization, such as the SR method and the steepest de- 
scent one, a basic ingredient is that at each iteration the 
new set of parameters {ce' k } are determined close enough 
to the previous set {ak} in terms of a prescribed dis- 
tance. The fundamental difference between the SR min- 
imization and the standard steepest descent method is 
simply related to the definition of this distance A a . 

Within the SR scheme, A Q is chosen to be the square 
distance between the two normalized wave functions cor- 
responding to the two different sets of variational param- 
eters {a' k } and i.e., 



A (SR) = 2 _ 2 



(*{a»}l*{a i} > 



(27) 



The reason to normalize the two wave functions before 
computing their distance is obvious because, within a 
VMC scheme, the normalization of the wave function is 



irrelevant for quantum mechanical averages. The basic 
advantage of the SR method is the possibility to work 
directly with the wave function distance Al SR '. In fact, 
this quantity can be explicitly written in terms of the ma- 
trix s [Eq, (25)] by substituting Eq. (15) in its definition 
Eq. (27), yielding' 

A i SR) = E E - «*)(«{ - «i) + 0(\a k a' k ?). 



fc=i i=i 



(28) 

Therefore the most convenient change of the variational 
parameters is to minimize the functional 

T({a k - a k }) = AE + AA^ SR) . 



Here AE is the linear change in the energy AE = 

~ Sfe/fe( Q! /s — an d f° r ^a R ' the leading term of 
the expansion in small a' k — a k given in Eq. (28) can 
be used. A is a Lagrange multiplier that allows for 
a stable minimization with small change Ai 31 ^ of the 
wave function [^{a*,})- Then the stationary condition 
5!F({a' k — ak})/S(a' k — a k ) = naturally lead to the SR 
iteration scheme described by Eq. (23) with St = 1/(2A). 

In a similar manner, it is also possible to obtain the 
standard steepest descent method. In this case the Carte- 
sian metric defined in the p-dimcnsional space of the vari- 
ational parameters is implicitly assumed to distinguish 
the two sets of variational parameters, i.e., the distance 
here is defined to be 



Ai SD > = ek - «*) S 



k=i 

The same argument used above to minimize AE + 
AA^a^ with respect to the variational parameter change 
ak — a' k will then lead to the standard steepest descent 
algorithm described by Eq. (17). 

The advantage of the SR method compared with the 
steepest descent one is now transparent. Sometimes a 
small change of the variational parameters corresponds 
to a large change of the wave function, and conversely a 
large change of the variational parameters can imply only 
a small change of the wave function. The SR method 
takes into account this effect through a better definition 
of the distance Ai SR) [Eq. (27)]. 
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Here a single SR iteration of the SR minimization 
scheme is summarized as follows: i) a set of variational 

parameters {ctk} — { a k } is given after the i-th iter- 
ation, ii) the generalized force fk [Eqs. (16) and (18)] 
and the matrix s [Eqs. (25) and (26)] are calculated 
statistically using a small variational Monte Carlo sim- 
ulation (bin) containing typically a few thousand sam- 
ples (M = 1000 4- 10000) distributed according to the 
wave function squared [^/^"[.(a;)! , and hi) a new set of 

the variational parameters {o:l <+1 ^} is determined from 
Eq. (23) with a suitable choice of St. After a few hundred 
(or sometimes thousand) iterations needed for equilibra- 
tion, the iteration described above is further repeated 
to statistically determine the optimized variational pa- 
rameters until the desired statistical accuracy is reached. 
This method allows for a very accurate determination 
of the optimized variational parameters with very small 
statistical uncertainty. In the present study, all the varia- 
tional parameters are optimized using this SR minimiza- 
tion method. 

A sample case study for the SR minimization scheme is 
presented in Fig. 2 for the ID antiferromagnetic Heisen- 
berg model on a L = 22 ring. As seen in Fig. 2, after 
the first few hundred iterations needed for equilibration, 
the variational parameters fluctuate around the stable 
mean values. It is interesting to notice in Fig. 2 that the 
variational parameters may continue to change substan- 
tially, even after the energy appears to reach its equi- 
librium value only after the first ~ 50 iterations. The 
reason is that a very tiny energy gain, not visible in this 
plot, is implicitly reached at equilibrium, by satisfying 
very accurately the Euler condition of minimum energy, 
i.e., fk = for k = 1,2, • •• ,p. A stochastic minimiza- 
tion method like the steepest descent method or the SR 
one, which is based not only on the energy but also on 
its derivatives (fk), is therefore much more efficient, es- 
pecially for statistical methods where energy differences 
for slight changes of the variational parameters are often 
very noisy. 

Having determined the optimized variational parame- 
ters {a* k } as described above, a single standard VMC sim- 
ulation is performed to calculate various physical quanti- 
ties for (^/{a*}) using a bin length much larger than the 
one used in the SR minimization procedure. Neverthe- 
less, the mentioned Euler conditions (fk — for all k) are 
usually satisfied within a few standard deviations simply 
because the statistically averaged variational parameters 
are much more accurate than the ones obtained at the 
final SR iteration. 26 



C. Effective Hamiltonian approach 

As is well known, within the variational approach, it is 
sometimes very difficult to describe long range properties 
accurately. This is simply because these long range prop- 
erties are rather insensitive to the energy, which is in- 
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FIG. 2: Monte Carlo evolution of (a) variational parameters 
(Ai and A3) and (b) energy as a function of Stochastic Recon- 
figuration (SR) iterations for the ID spin-1/2 antiferromag- 
netic Heisenberg model Hid on a L — 22 ring. Here the SR 
method with St = 0.2/ J [Eq. (23)] is used. In (a), Ai and A 3 
are the two variational parameters of the wave function (for 
details, see in IIIB), which are both initialized to the value 
Ai = A3 = 1. In (b), at each iteration, the energy is com- 
puted by a small VMC simulation for the wave function with 
the fixed variational parameters given at the corresponding 
iteration in (a). 



stead determined mainly by short range correlation func- 
tions. A clear example of this limitation is found when 
this method is applied, for instance, to study the long 
range antiferromagnetic order for the spin-1/2 antiferro- 
magnetic Heisenberg model on the square lattice. The 
long range antiferromagnetic order is now considered a 
well accepted and established property of this model. 27-29 
Nevertheless it has been shown in Ref. 30 that two wave 
functions with completely different long-distance proper- 
ties, with or without antiferromagnetic long range order, 
provide similar and very accurate (within 0.1% accuracy) 
energy per site in the thermodynamic limit. 

Here we shall consider a possible way to overcome this 
limitation of the variational method by introducing what 
will be named "effective Hamiltonian approach" in the 
following. The main idea of this method is to approx- 
imate, as accurately as possible, the exact Hamiltonian 
H by an effective Hamiltonian H cS , for which the ex- 
act ground state can be numerically sampled by Green 
function quantum Monte Carlo schemes. The important 
point is that, within this construction, some of the im- 
portant short range properties of the Hamiltonian H are 
preserved, and therefore a much better control of correla- 
tion functions is achieved compared to a mere variational 
approach. 

The short range properties of the Hamiltonian H can 
be conveniently defined in the configuration basis {x}. 
Each element of this basis is defined by the positions 
and spins of the N electrons. In this basis the ma- 
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trix elements of the Hamiltonian H are then indicated 
by (x\H\x') = H x x i. In this paper, we define a short 
range Hamiltonian as a Hamiltonian H for which the 
off-diagonal matrix elements H xx i are non zero only for 
configurations x and x 1 connected one to the other by 
local short-range "moves" of electrons (e.g, hoppings or 
spin- flips). More precisely, 



^ for \x-x'\ < R 
= otherwise 



(29) 



where R determines the range of the off-diagonal matrix 
elements. In this definition, not only conventional Hub- 
bard, Heisenberg, and t-J models can be thought of short 
range Hamiltonians (with R = 1), but also models with 
any interaction V which is diagonal in the configuration 
basis {x} {V XyX i cx V X S X:X '), e.g., the long range Coulomb 
interactions, can be considered short range Hamiltonians. 
Therefore, within this definition, most of the physically 
relevant Hamiltonians arc short range, essentially be- 
cause quantum fluctuations are important only through 
the short range off-diagonal matrix elements, in the ab- 
sence of which the Hamiltonian is simply classical. 

At present, short-range spin Hamiltonians which can 
be solved numerically by quantum Monte Carlo methods 
are the ones for which the quantity 



s x ', x = ip G {x')H x , !X ip G (x) < 



(30) 



is strictly negative for all non zero matrix elements 
H X '. X . Whenever a wave function iPg(x) satisfying 
this condition exists, a unitary transformation H x x > — > 
Sgn [^g(x)] Sgn [ipoix')] H XjX i transforms the Hamilto- 
nian to a standard bosonic one (namely with all non- 
positive off-diagonal matrix-elements) that can be solved 
numerically with quantum Monte Carlo techniques, with- 
out facing the so-called "sign problem". 31 For instance, 
for the antiferromagnetic Heisenberg model in any bi- 
partite lattice, a particularly simple variational wave 
function {x\4>g) = 4>g{x) satisfying the Marshall sign 
rule, tpo(x) cx (— 1)^*, allows to satisfy the condition 
s x .x' < 0, where Na is the number of electrons with 
down spin on one of the sublattices for the given configu- 
ration {x} (see App. B). With this wave function 
it is therefore possible to solve numerically the spin- 1/2 
antiferromagnetic Heisenberg model in ID, on a two-leg 
ladder, and on the 2D square lattice. In these cases, it is 
also clear that with the same Marshall sign of the wave 
function ipa(x), different low energy properties can be 
obtained, i.e., a gapless spin liquid ground state for the 
single chain, a gapped spin liquid for the two-leg ladder, 
and a quantum antiferromagnet for the 2D limit. 

However, as is well known, only for very particular 
models the Marshall sign rule and Eq. (30) are satis- 
fied, 28,29 ' 32 and in general Eq. (30) is violated regardless 
of the wave function \4>g) used. Even when a wave func- 
tion \iPg) with the optimal signs, i.e., the exact GS signs, 
is used for generic frustrated Hamiltonians, there still ex- 
ist off-diagonal matrix elements with s XtX > > (notorious 



sign problem). Namely, due to frustration, they do not 
decrease the energy expectation value. 

In order to overcome this difficulty and treat more 
generic models, we will define below an effective Hamil- 
tonian H cS , which is closely related to H, by using an 
optimal wave function ipa (%) , that we name in the follow- 
ing the " guiding function" . This guiding wave function is 
required to be non zero for all configurations {x}. Even if 
this requirement is not satisfied, all the forthcoming anal- 
ysis remain valid in the subspace of configurations {x} for 
which ipa{ x ) 7^ 0. Once the guiding wave function is pro- 
vided and thus H cS is defined, the effective model system 
H cS is solved exactly using the standard Green function 
quantum Monte Carlo method. 33 As will be shown later, 
the low energy properties of H are weakly dependent 
on the low energy properties, i.e., long range behavior, 
of ipc(x). 

The GS wave function of H is approximated by the 
GS of an effective Hamiltonian H cS . This approximate 
variational state is very good in energy because most of 
the matrix elements of H are treated exactly in H cil , 
whereas the remaining ones are removed and traced to 
the diagonal terms of H cR . As it will be shown later 
on, this enforces the constraint that the GS of H eS has 
the same non trivial signs of ^Pg( x )- Therefore, if ipG{x) 
is chosen to have the same signs of the exact ground 
state of H for most configurations, this approach becomes 
essentially exact. 34 

Within the effective Hamiltonian approach, the prob- 
lem to construct an accurate approximate wave function 
for the GS of H is therefore related to how to know phases 
of the GS. This, we believe, is a much simpler task, be- 
cause a good variational wave function of a short range 
Hamiltonian should provide also good phases. In fact, 
the variational energy of the guiding function |V>g), 



E{^ g ) 



T, x , x > s S n [sx,x'] \il> G (x)ip G (x')Hx, x 



depends strongly on the signs of il> G {x) via the short 
range terms appearing in s x , x '- This assumption that 
a variational wave function with accurate energy should 
have very good signs can be checked directly on small 
size clusters. Moreover, for the nearest neighbor anti- 
ferromagnetic Heisenberg model on the square lattice, it 
is known that all good variational wave functions satisfy 
the Marshall sign rule, although they may display differ- 
ent large distance behaviors. 30 This fact clearly confirms 
this assumption because all these wave functions differ 
only in their amplitudes but not in their phases. 

From these considerations, in this paper, we will chose 
as a guiding function the projected BCS wave function 
|p-BCS) described in the previous subsection with the 
variational parameters optimized for each J' / J by a care- 
ful energy minimization using the SR scheme (Sec. II B). 
In most cases, due to the quality of the guiding function 
used, no relevant corrections are found in the various 
large distance correlations calculated for |p-BCS) when 
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these are compared with the ground state correlations 
of H eS . In some special cases, the use of the effective 
Hamiltonian approach is instead of crucial importance. 



1. Definition of the effective Hamiltonian 

As mentioned previously, the effective Hamiltonian 
H cS is defined in terms of the matrix elements of H, 
which are chosen to generate a dynamic as close as pos- 
sible to the exact one. An obvious condition to require is 
that if \ip G ) is exact, then the ground state of H eS and its 
eigenvalue have to coincide with the ones of H. In order 
to fulfill this condition, the so-called lattice fixed node 
(FN) was proposed 35 . In the standard lattice fixed node 
approach, all the matrix elements which satisfy Eq. (30) 
are unchanged, whereas the remaining off-diagonal ma- 
trix elements are dealt semiclassically and traced to the 
diagonal term, defining the standard FN Hamiltonian 
jjes = ^fn Thc FN Hamiltonian £fn ig obtained by 

modifying its diagonal elements in a way that the local 
energies corresponding to the FN Hamiltonian H FN and 
the exact one H coincide for all configuration x, namely, 



(jPg\H\x) = (iPg\H \x) 
(iPg\x) {^g\x) 

and therefore the FN Hamiltonian is defined by 
H 



(31) 



H x\x - t 



where 



, , if x' 7^ x and s x ^ x < 

if x' 7^ x and s x >, x > 

H x . x + V st {x) i{x' = x 

(32) 



V sf (x) = ]T yj G (x')H xl . x /yj G (x). 

{x'{^x),s x , iX >0} 

The FN approach was inspired from the similar fixed 
node method on continuous systems, 16 and indeed is a 
well established approach 35 which provides also varia- 
tional upper bounds of the ground state energy, i.e., 
E$ N < E(ip G ) where is the ground state energy 

of H™. 

As we will show below, for lattice systems there is a 
better way to choose this effective Hamiltonian 36 , which 
not only provides better variational energies, but also al- 
lows for a better accuracy of the low energy long distance 
properties of the ground state. For this end, it is impor- 
tant to notice the following key difference between a lat- 
tice system and a continuous one: in the lattice system 
the configurations {x} which do not satisfy the condition 
(30) may be a relevant fraction of the total number of 
configurations, whereas in the continuous case such con- 
figurations represent just an irrelevant "nodal surface" of 
the phase space. Therefore dropping all the off-diagonal 
matrix elements with s XtX i > as in the standard FN 
method seems to be a too drastic approximation for the 



lattice case. This approximation can be indeed improved 
for lattice systems because, contrary to the continuous 
case, the FN Hamiltonian does not provide the best vari- 
ational state with the same signs of the chosen guiding 
function ipc(x)- 

The main consequence of neglecting all the off-diagonal 
matrix elements with s XjX i > is a bias of the dynamic, 
as electrons cannot move freely in some of the configura- 
tions. In order to compensate this bias in the diffusion 
of the electrons, we introduce a renormalization constant 
K < 1, which reduces the off-diagonal "hopping" in the 
allowed configurations with s x , x > < 0: 



l' ,X * 



KH X > :X if x' 7^ x and s x >, x < 

if x' x and s x >. x > 

H x . x + V sf (x) Xxf = x 

(33) 

whereas in order to satisfy the condition (31) V sf is mod- 
ified as follows: 

V si (x) = (1-K) Yl ^ G (x')H x , tX /ip G (x) 

{x'(^x),s x , >x <0} 

+ J2 ip G {x')H x , !X /ip G (x). 

{x'&x),s x , <x >0} 

2. Optimal choice for the constant K 

Whenever there is no sign problem and s XyX > < 0, 
K = 1 is obviously the best choice for which H = H cR . 
Also this choice K = 1 coincides with the standard lat- 
tice FN approach [Eq. (32)] (which will be denoted by 
thc acronymous FN). Conversely, when the sign of thc 
off-diagonal term in H are frustrated {s XjX > > 0), a better 
choice of the constant K can be obtained by using a rela- 
tion which has been well known for continuous systems, 
and used to correct efficiently the error due to the fi- 
nite time slice discretization in the diffusion Monte Carlo 
(DMC) calculations. 37 

Let us first discuss this relation used in DMC before 
considering the lattice case. In the DMC, the small 
imaginary time (At) evolution of the electron configu- 
rations is governed by the exact Hamiltonian propaga- 
tion, \ip ) — > cxp(— HAr )\yj), with a diffusion coefficient 
D determined only by the free kinetic operator in H. It 
is then possible to correct the approximate finite At dy- 
namic corresponding to the fixed node Hamiltonian, by 
requiring that it satisfies exactly the short time diffusion 
condition. This condition can be simply written as 



[x, [H,x]] = D, 



(34) 



where D = 3h jm is the diffusion coefficient in three 
dimensions, and x is the electron position operator. 

For lattice systems with periodic boundary conditions 
(PBC), the position operator x is not well defined, as it 
cannot be matched with the boundary conditions. This 
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limitation can be easily solved by using periodic spin po- 
sition operators defined in the exponential form by 



exp 



R 



{hu ■ R) S% 



(35) 



where S z s is the ^-component of the spin operator at 

site R, and v labels the spatial coordinates, e.g., h x = 
(2ir/l, 0) and h y = (0, 2n/l) for a L — I x I square lattice. 
These operators are diagonal in the basis of configura- 
tions {x}, {x\X v \x') = X„(x)5 x . x >, as is the analogous 
position operator x in the continuous case. Remarkably 
X v is exactly equivalent to the well known Lieb-Schultz- 
Mattis operator, used to show a well known property on 



J 



the low energy spectrum of spin-1/2 Heisenberg Hamilto- 
nians 38 . After simple inspection, a relation similar to the 
one (34) can be found also for the periodic spin position 
operators X v , by simply imposing the following equation: 



X 



H, X v 



Hg) = (ipc 



X 



H cjs ,X v 



VI>g) 

(36) 

For the lattice case, both the left hand side and the right 
hand side of this relation have non trivial expectation val- 
ues. Both of them can be simply calculated by a standard 
VMC method without particular difficulty, after recast- 
ing them in a more conventional form for VMC calcula- 
tions: 



xl 



H, X v 



= -J2J2^G(x)^ G (x')n x , xl \Xu(x)-X u (x')\' 



x x > 



Y J ^G{x')n XtXl \X v {x)-X u {x l )\ 2 /^ G {x) , (37) 

. x' 

I 



a relation that is obviously valid both for H = H and 
for TL = H cS . Here we assumed W XiX / = W x ',x- There- 
fore, similarly to the continuous scheme 37 , the value of 
the constant K can be determined from Eq. (36) with 
high statistical accuracy. Notice that if there is no sign 
problem, i.e., s x , x ' < for all configurations {x}, the 
constant K turns out to be exactly one with vanishing 
statistical error, yielding again H cS = H. 

After determining the constant K, the effective Hamil- 
tonian H cS [Eq. (33)] is defined, and the ground state 
\ifjo S ) with its eigenvalue Eq S and the corresponding low 
energy excitations of H cS can be computed using the 
standard Green function quantum Monte Carlo method 33 
without sign problem. 

To compute the expectation value of the energy 



(% S \H\% S ) 



(38) 



using \ipQ ) as an approximate ground state for H, the 
method described in Rcf. 39 can be applied, which very 
often improves sizably the upper bound estimate of the 
energy, i.e., E(tpQ S ) < Eq S < E(tpc), even in the stan- 
dard FN case with K = l. 37 As also remarked in Ref. 39, 
contrary to the continuous case, for lattice Hamiltonians 
the lowest variational energy E(tpQ S ) does not correspond 
to K = 1 in general. 

In the following, we will indicate by FNE the improved 
FN effective Hamiltonian method [Eq. (33)] with the con- 
stant K < 1 determined by the condition (36), using for 
|V>g) the lowest energy variational wave function of the 



form described in Sec. II A. 



D. Calculation of dynamical correlation functions 

The spin-one excitation spectrum of the effective 
Hamiltonian H cS can be calculated by applying the for- 
ward walking technique 33 used to evaluate the imaginary 
time evolution of the following quantity 



S(k,r) 



-tH" 



SIM") 



{ip G \e- 



rff'ff |,/,efl: 



(39) 



where = ^= ^ r e tk r S* and IV'o 9 ) is the ground state 

of H cS . Note that the imaginary time propagation in 
Eq. (39) can be evaluated without discretization errors 
in time r, as is common to many other Quantum Monte 
Carlo techniques, 40 and also pointed out in Ref. 41 for 
the present Monte Carlo scheme. By simple inspection, 
the spin one excitation energy E^ =1 of H cS for mo- 
mentum k can be calculated by fitting the large imag- 
inary time behavior of S(k.,r) oc exp [— _E(k)r] . Here 
E(k) = E^ =1 - E^ B and E"f is the ground state energy 
of H cB . Thus we fit logS(k,r) with 



log 5(k, r) = — r£(k) +A + B log(r) 



(40) 



for t > r c , where A, B, and E(k) are fitting parameters, 
and r c is a suitable cutoff time, large enough so that the 
fitting form (40) can be used with good accuracy. In 
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fact, the present fit is very stable in r, and a satisfactory 
convergence of E(k) is obtained even for relatively small 

T C ~ 2/ J. 

As a typical example, Fig. 3 (a) shows a semi log plot 
of S(k,r) at k = w for the ID spin-1/2 antiferromag- 
netic Heisenberg model Hid- It is clearly seen that for 
a wide region of r > 2/J, log S(k,r) is linear so that 
we can safely estimate the excitation energy E(k) with 
Eq. (40). The error bars for E(k) can be efficiently eval- 
uated by using the "bootstrap technique" , 42 because the 
numerator and the denominator of Eq. (39) are highly 
correlated. To demonstrate the accuracy and reliability 
of the method, the calculated E(k) for Hid is shown in 
Fig. 3 (b). As seen in the figure, the agreement between 
our results and the exact values is excellent. 

It is worthwhile to emphasize that, within this method, 
a single Monte Carlo simulation allows to calculate all the 
lowest spin one excitation energies for various momenta 
with no extra effort. Moreover, the method is not re- 
stricted to the computation of the spin triplet spectrum 
alone, but can be easily generalized to arbitrary operators 
O as long as they are diagonal in terms of the chosen basis 
{x}. Although the high energy excitations are more dif- 
ficult to calculate because the signal decays much faster 
in time, the most important low energy spectrum can be 
accurately determined. 




FIG. 3: (a): Imaginary time r evolution of S(k,r) at k = ir 
for the ID spin-1/2 antiferromagnetic Heisenberg model Hm 
on a L = 22 ring, (b) Excitation energy E(k) = E% =1 - E<f 
extracted from S(k,r) with large r shown in (a) (squares). 
The error bars are smaller than the size of the symbols. For 
comparison, the exact values are also presented by solid cir- 
cles. 



III. NUMERICAL RESULTS 

A. Projected BCS wave functions and particle-hole 
symmetry 

Before presenting our numerical results, here we will 
show an important relation between the projected BCS 



wave function [p-BCS) and the particle-hole symmetry, 43 
which turns out to be crucial to differentiate spin systems 
defined on the conventional non-frustrated square lattice 
and the ones defined on the triangular lattice geometry. 
In both cases each site of the lattice is represented by 
a vector r = T\Ty + t^t-i with v\ and ri being integer 
(see, e.g., Fig. 1). Then a particle-hole transformation is 
defined with 

4, CT - (-l) ri+r2 Cr, CT , (41) 

or in the reciprocal lattice space with the reciprocal lat- 
tice vectors g*\ and 52 , this transformation is equivalently 
defined with 

C L -» c -k+Q, CT (42) 

where Q = (g± + (fe)/2. As shown in App. B, when- 
ever the BCS Hamiltonian Hbcs is invariant under the 
particle- hole transformation and A& is real, i.e., 

!£k = -£-k+Q 
A k = -A_ k+ Q (43) 
M = 0, 

the corresponding projected BCS wave function |p-BCS) 
[Eq. (8)] satisfies the so-called Marshall sign rule. 44 The 
Marshall sign rule is an exact property for the ground 
state of the spin-1/2 antifferomagnetic Heisenberg mod- 
els on non frustrated lattices such as the square lattice, 
and indeed for these models the minimum variational en- 
ergy is achieved when this rule is satisfied by the pro- 
jected BCS wave function |p-BCS). 

B. One dimensional limit and spin fractionalization 

We shall first show the results for the uncoupled chain 
limit, i.e., for the ID spin-1/2 antiferromagentic Heisen- 
berg model with nearest neighbor coupling: 

H 1D = J^Si-Sj. (44) 

(id) 

Several previous studies 18,45 have found that the ground 
state of Hid can be described very accurately by the pro- 
jected BCS wave function |p-BCS) [Eq. (8)] with only 
the nearest neighbor hopping ti,j = #ij±i and the first 
three neighbors for the gap function Aj„- = Aidij±i 
(I = 1,2,3). Here the site-i is represented by the vector 
ri = i T\ . Notice that the ground state of H\d satis- 
fies the Marshall sign rule, and therefore from Eq. (43) 
with Q = 7T, [i and A; with I even have to be identically 
zero. It was also pointed out 18 that the inclusion of the 
third neighbor gap function A3 is crucial to improve the 
accuracy of the |p-BCS). As shown in Fig. 2(a), the op- 
timized parameters for L = 22 are Ai = 2.947 ± 0.003 
and A 3 = 0.737 ± 0.002 with the notations of Eq. (7), 
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i.e., Afe = 2Ai cos fc + 2A3 cos 3fc. The correspond- 
ing variational estimate of the total energy is E/J — 
—9.78411 ± 0.00005 which is in excellent agreement with 
the exact value of E/J — —9.78688. For larger clusters 
the variational parameters change slightly and smoothly, 
and the same kind of accuracy is obtained even in the 
infinite size limit. By simple quadratic (linear) extrapo- 
lation in 1/1/ for the energy (A;) with data up to L = 150 
sites, we found E/JL = -0.442991(3) [Ai = 3.41(3), 
A 3 = 0.90(1)] which compares very well with the well- 
known exact value E/JL = -(In 2 - 1/4) = -0.443147. 
This result suggests that our variational method is par- 
ticularly accurate in describing this ID exactly solvable 
case, which is precisely our starting point for studying 
weakly coupled chains with J' =/= 0. 

In the rest of this subsection, we shall show that also 
the low-energy excited states can be described by pro- 
jected BCS states. As we mentioned above, to satisfy 
the Marhsall sign rule, the optimized variational param- 
eters in |p-BCS) for the ground state meet the constraint 
relation (43) with Q — it. Therefore, at this momen- 
tum k = Q, the spin-1/2 BCS excitation spectrum E k 
[Eq. (10)] shows gapless excitations at fc = ±Q/2, in 
perfect agreement with the exact spinon spectrum of the 
Bethe-ansatz solution. 14 Since the elementary excitations 
of the BCS Hamiltonian with energy E^ are described by 
the standard Bogoliubov modes, 46 



7-fc,J, = w fc c l,T + u kC-k,i, 



(45) 



the simplest variational state for the spinon at momen- 
tum k is 



|k)=7>G7ljBCS) 



(46) 



To see whether this state |k) corresponds to a spinon 
state, we consider a ring with odd number of sites L = 31 
and z-component of the total spin S* ot — —1/2. For this 
case it is known that a well defined spinon exists only for 
half of the total Brillouin zone (tt/2 < |k| < tt). 47 For 
this branch, as shown in Fig. 4, the wave function |k) 
represents very well a spinon with momentum k, as can 
be verified by the good accuracy in the energy and its 
small variance 



^ 2 (|k)) = 



(k|g 2 |k) 
(k|k) 



(k|g|k) 
(klk) 



(47) 



The variance <7 2 (|k)) is zero for an exact eigenstate, and 
is small for a very accurate variational state. As seen in 
Fig. 4, this is clearly the case for the momenta fc > tt/2. 

For the remaining branch of momenta (|k| < tt/2), the 
excitations for Hid are no longer elementary. 47 As shown 
in Fig. 4, although the state |k) is formally defined even 
for those momenta, it represents a poor representation for 
the exact excited state. In fact -E'(k) becomes quite off 
from the exact value, and a 1 (|k)) considerably increases 
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FIG. 4: Lowest spin-1/2 excitation energy E(k) and E(k ( - 3 ') 
(lower panel) and the variance a 2 (\k)) and a 2 {\k^)) (upper 
panel) of the projected BCS states with a single BCS mode 
jfc) and three BCS modes \k^) (see text) as a function of 
momentum k. The model studied is the ID spin-1/2 anti- 
ferromagnetic Heisenberg model on a ring with L=31. For 
comparison, in the lower panel, the exact lowest excitation 
energy is plotted by crosses, and the BCS spectrum de- 
noted by dashed line is scaled by a factor to match the exact 
bandwidth. The remaining lines are guides to the eye. 



when the momentum crosses the " spinon Fermi surface" . 
Remarkably, by projecting a non-elementary excitation 
of the BCS Hamiltonian, namely, 



|k (3) )=^G7l a 7l kFa 7l, T |BCS) 



(48) 



with kp = tt(L + 1)/2L, which mimics the correct 3- 
spinon eigenstate, we can achieve a good agreement for 
the spectrum even in the region outside the spinon Fermi 
surface (see Fig. 4). Notice that, although the projection 
Vq is crucial to gain a quantitative agreement for the 
spectrum, the BCS spectrum E^ already gives a quali- 
tatively correct feature of gapless excitations with finite 
spinon velocity at the right momentum k = tt/2 (see 
Fig. 4). 

Since the state |k) [Eq. (46)] is very accurate only in 
the momentum region 7r/2 < |fc| < 7r, we conclude that 
only half of the elementary excitations of the BCS Hamil- 
tonian remain to be well defined excitations for the spin 
Hamiltonian Hid after applying the Gutzwiller projec- 
tion Vq. This effect is expected to hold also in a 2D 
fractionalized phase. It turns out that the elementary 
BCS excitations which describe the correct spinons after 
applying the projection operator Vq can be obtained by 
adiabatically switching off the gap function Afe, a process 
that defines naturally a Fermi surface, so that fermion 
quasiparticles can be created (destroyed) in unoccupied 
(occupied) states only. 
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In 2D it is very difficult to confirm directly the above 
"selection rules" of the Gutzwiller projection acting on 
the elementary BCS excitations. For instance, to study 
single spinon excitations, one might think of a 2D sys- 
tem on (I x I) lattice with I odd as a natural extention of 
the ID system considered above. However this 2D system 
should contain many elementary excitations of spinons as 
it is easily understood by considering the J' — > limit, 
and thus it is not an ideal system for studying single 
spinon excitations. Nevertheless, it is very important to 
have reached a very accurate ID limit, with the correct 
spin fractionalization, within the present variational ap- 
proach, because this approach can be easily extended to 
higher dimensions. 

Finally, it is worth mentioning that, according to the 
gauge theory by Wen, 48,49 the low energy excitations of 
several 2D spin liquids described by this |p-BCS) varia- 
tional ansatz, can be understood at the mean field level 
without taking into account the Gutzwiller projection, 
simply because this projection becomes irrelevant for 
large distance correlations. This is the case for a "Z 2 
gapless spin liquid". Without entering into too much 
details of this theory, we only mention here that a Z 2 
spin liquid can be described by a |p-BCS) for which 
no SU(2) gauge transformation — remaining in the re- 
stricted Hilbert space with no doubly occupation — al- 
lows to eliminate the gap function in the corresponding 
BCS Hamiltonian. In such a case, only the gauge trans- 
formation Cj j(7 — > —c i a leaves the BCS Hamiltonian un- 
changed, a transformation which defines the Z 2 group 
toghcthcr with the identity operation. Most of the spin 
liquids which we will describe in this paper are Z 2 spin 
liquids in the triangular lattice geometry and should be 
stable according to the theory mentioned above 48,49 



C. Weakly coupled chains in the triangular 
geometry with J'/J small 

In this subsection, as a first step towards the 2D limit, 
we shall consider the spin-1/2 antiferromagnetic Heisen- 
berg model for weakly coupled chains in the triangu- 
lar lattice geometry with J'/J small (see Fig. 1). Note 
that this region is appropriate for the material Cs 2 CuCLt, 
where J'/J ~ 1/3. 11 



1. Variational results 

Motivated by the great success of the present varia- 
tional approach in the ID system discussed in the pre- 
ceding section, the variational ansatz state which we shall 
consider here for this 2D model in the region J'/J small 
is a similar projected BCS state |p-BCS) described by 
Eqs. (4)-(10). Here both t itj and A itj in the BCS Hamil- 
tonian [Eq. (4)] parameterizing |p-BCS) are assumed 
translational invariant, and therefore they depend only 
on the relative vector I = fi — fj = l\f\ + I2T2 = {h, h) 



(li and I2: integer) between the two sites i and j, i.e., 
and Ap Also the C 2v point group symmetry is as- 
sumed for the variational parameters, e.g., A ; -= A K ; - = 
A K j*= A K K j- for the A\ symmetry, where TZ X (1Z V ) is 
the reflection operator about the xz (yz) plane: 

TZJ = (h + l 2 )n - 1 2 t 2 
Hyl = -(h + k)Ti + Z2T2. 

To optimize the variational wave function, we use the SR 
optimization method described in Sec. II B, which en- 
ables us to determine the optimized variational parame- 
ters with very high accuracy. 

As shown in Tab. I, the optimized tj- is found to be 
non zero only for the nearest neighbors along the chain 
direction fi (wc set £5 = 1), and negligible otherwise, 
whereas the optimized A ; - is instead found to be sizable 
even among different chains (for instance, A/ 2l i), shown 
in Tab. I), displaying a true two dimensional character. 
It is also found that the symmetry of the gap function A,- 
which minimizes the variational energy is the A\ symme- 
try under the C 2v point group. Finite-size corrections to 
the variational parameters scale as 1/L, and the 18 x 18 
cluster is found large enough to be close to the thermo- 
dynamic limit, at least, for not too small values of J'/J. 
All the variational parameters for the 18 x 18 cluster as 
a function of J'/J(< 0.7) are tabulated in Tab. I along 
with the variational energy. It should be noted that since 
the Marshall sign rule is no longer satisfied for this model, 
the constraint relation (43) does not have to be met by 
the variational wave function. For example, as seen in 
Tab. I, the optimized chemical potential n turned out to 
be different from zero. 

To gain better insight on the physical nature and prop- 
erties of these variational states, let us next evaluate the 
BCS excitation spectrum given by Eq. (10) in the 
thermodynamic limit with the optimized variational pa- 
rameters. In Fig. 5 and Fig. 6, the contour lines for 
£k = ek — n = and Ak = are plotted together with 
the boundary of the first Brillouin zone (BZ) of the trian- 
gular lattice for J'/J = 0.33 and 0.5, respectively. 50 As 
can be seen clearly from these figures, £k and Ak van- 
ish at the same momenta with incommensurate "nodal" 
points, and t hus the co rresponding BCS excitation spec- 
trum £"k = \/£k + snows gapless excitations at these 
momenta. It is also interesting to notice that with in- 
creasing J'/J the shape of the contour line of the gap 
function Ak = changes form open lines to a closed one 
in the BZ, acquiring a clear two dimensional characteris- 
tic. On the contrary, the minimum variational energy is 
always achieved with negligible values of the inter chain 
hopping t^even for the largest J'/J considered. 

Before considering the 2D limit, it is important to dis- 
cuss in more details the properties of few coupled chains 
in the triangular lattice geometry displayed in Fig. 1. If 
we consider only two chains in the fi direction, this sys- 
tem corresponds to the well known zig-zag ladder with 
couplings J and J'. It is now well established that 



13 



3'/3 


0.1 


0.2 


0.33 


0.4 


0.5 


0.6 


0.7 




-0.066(4) 


-0.085(5) 


-0.109(4) 


-0.147(4) 


-0.200(4) 


-0.232(4) 


-0.253(3) 


A(i, ) 


2.28(2) 


2.02(2) 


1.74(2) 


1.61(3) 


1.48(3) 


1.36(2) 


1.24(3) 


A( ,i) 


0.545(7) 


0.617(6) 


0.627(6) 


0.644(8) 


0.689(8) 


0.739(7) 


0.787(7) 


A d,i) 


0.155(4) 


0.179(3) 


0.192(3) 


0.214(4) 


0.256(4) 


0.293(3) 


0.329(4) 


A(„i, 2 ) 


0.007(5) 


0.099(6) 


0.147(6) 


0.160(6) 


0.162(8) 


0.153(4) 


0.143(6) 


A(2,0) 


-0.010(1) 


-0.006(2) 


-0.0012(8) 


0.002(2) 


0.013(1) 


0.025(2) 


0.037(2) 


A(0,2) 


0.195(3) 


0.159(5) 


0.070(3) 


0.058(4) 


0.062(4) 


0.068(5) 


0.066(4) 


A( 2 ,l) 


0.346(4) 


0.381(4) 


0.370(4) 


0.363(5) 


0.360(5) 


0.361(6) 


0.356(6) 


A(l,2) 


0.041(2) 


0.085(3) 


0.098(3) 


0.102(4) 


0.102(5) 


0.097(6) 


0.095(4) 


A(-2,3) 


0.052(1) 


0.040(1) 


0.006(2) 


0.006(2) 


0.010(2) 


0.013(2) 


0.015(1) 


01 


0.491(5) 


0.434(4) 


0.383(4) 


0.363(5) 


0.345(5) 


0.327(6) 


0.304(6) 




1.00243(3) 


1.01014(7) 


1.0286(2) 


1.0434(2) 


1.0700(3) 


1.1030(3) 


1.1425(3) 


-Evmc / 3L 


-0.44590(1) 


-0.44687(1) 


-0.44929(2) 


-0.45118(2) 


-0.45474(2) 


-0.45932(2) 


-0.46514(2) 


Efn/ 3L 


-0.446074(1) 


-0.44723(1) 


-0.45005(1) 


-0.45223(1) 


-0.45628(1) 


-0.46156(1) 


-0.46823(1) 


-Efne 1 3L 


-0.446075(1) 


-0.447233(2) 


-0.45007(1) 


-0.45229(1) 


-0.45642(1) 


-0.46183(1) 


-0.46875(1) 



TABLE I: Optimized variational parameters of the wave function |p-BCS) for the spin-1/2 antiferromagnetic Heisenberg 
model on the anisotropic triangular lattice [Eq. (1)] with various 3' / 3 and L — 18 x 18. A( n m ) is the gap function A? for 
f = nfi +mT2- All other variational parameters are zero except for the chemical potential /i and the nearest neighbor hopping 
in the f\ direction, which is chosen to be one. The value of K, variational energy -Evmc = E(p -BCS), FN (FNE) ground state 
energy Ef~n ~ E™ (Ufne = Eq) with \ipa) = |p - BCS) are also presented. The number in parentheses is the statistical error 
bar of the quantity corresponding to the last digit of the figure. 
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FIG. 5: Loci of k-points where £k = (dashed lines) and 
Ak = (solid lines) determined by using the optimized varia- 
tional parameters for 3' / 3 = 0.33. The boundary of the first 
Brillouin zone (BZ) for the triangular lattice is also denoted 
by long dashed lines. 



the ground state of the zig-zag ladder is spontaneously 
dimerized with dimers along the rungs (ji direction) for 
J' I J < 4, and that the low-lying spin excitations are 
gapped. 51 For small values of J' /J, it is impossible at 
present to detect the dimer order parameter numerically 
on finite size clusters simply because this quantity van- 
ishes exponentially for J' / 3 — > 0. However, it was also 
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FIG. 6: The same as in Fig. 5 but for J'/J = 0.5. 



shown previously 52 that, for systems with finite num- 
ber of chains, projected BCS wave functions can show 
spontaneous dimerization provided there is a gap in the 
corresponding BCS excitation spectrum E^. This is man- 
ifestly the case for the zig-zag ladder system, as well 
as for any system with finite even number of chains, as 
the corresponding finite discretization for the momenta 
in the y-direction do not allow to fulfill simultaneously 
the two conditions Ak = and £k = and thus the 
presence of a gap in the BCS excitation spectrum is im- 
plied. Therefore, the projected BCS wave functions can 
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naturally describe the crossover from a finite number of 
coupled chains, dimerized and gapped, to a gapless and 
fractionalized spin liquid in 2D, implying that the limit 
of infinite number of chains may be highly non trivial for 
a spin liquid wave function. 

In order to understand the properties of the present 
projected BCS wave function |p-BGS), we now report 
several physical quantities. Fig. 7 shows the spin-spin 
correlation functions 



C{t) = 



(V#> 



(49) 



with \ip) = |p-BCS) calculated in the ti and T2 directions 
for typical couplings with J'/J = 0.33 and J'/J = 0.5. 
As seen in Fig. 7, even though the inter-chain spin cor- 
relations are very weak, the intra-chain spin correlations 
are appreciably large, at least, for the clusters studied. 
In order to examine whether magnetic long range order 
occurs, we have studied the system size dependence of 
C(/ max ri) (— P s ) at the maximum distance (/ max ) in the 
f\ direction compatible with periodic boundary condi- 
tions. There exists long range magnetic order only when 
P s is finite in the thermodynamic limit. Our results of 
P s are shown in Fig. 8 (a) for J' / J = 0.33 and clusters 
up to L = 42 x 42. From the finite size scaling presented 
in Fig. 8 (a), it is clearly concluded that the projected 
BCS state |p-BCS) is not magnetically ordered. This 
is apparently expected because the projected BCS wave 
function |p-BCS) is a spin liquid state in 2D. 

Another important quantity to study is the dimer- 
dimer correlation functions along the chain direction: 



D(l) 



M ( 5555, ¥ ) (s z _ r _ 



(50) 

As is well known, D(l) shows undamped oscillations at 
large distance \l\ — > oo when there is a spontaneous bro- 
ken translation symmetry characterized by a dimerized 
spin-Pierls phase. Since for systems with finite num- 
ber of chains projected BCS wave functions show dimer 
order, 52 it is not surprising to sec large oscillations in 
D(l) as a function of the distance, as shown in Fig. 9 for 
J'/ J = 0.33. In order to rule out the possibility of a finite 
dimer order for the present projected BCS state in 2D, 
the system size dependence of the dimer order parameter 
Pd is studied in Fig. 8 (b) for clusters up to L = 42 x 42. 
The square of this order parameter (Pj) is related to the 
calculated largest distance dimer correlations through 



Pj = 18 



[D(l/2n)-D((l/2-l)n)} 



for a cluster of L = I x I sites. 53 As seen in Fig. 8 (b), it is 
clear that Pd — » as L — > oo . As expected, this projected 
BCS wave function |p-BCS) does not show spontaneous 
dimerization in 2D, and therefore it represents a genuine 
spin-liquid wave function. 
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FIG. 7: Spin-spin correlation functions C(f) with r = hri + 
1 2 t 2 for 3' jj = 0.33 [(a) and (b)] and J'/J = 0.5 [(c) and 
(d)] calculated using both variational Monte Carlo (VMC) 
and effective Hamiltonian (FNE) methods. The cluster sizes 
used are 30 x 30 [(a) and (b)] and 18 x 18 [(c) and (d)]. Ci(l) 
and C2(l) indicate C(iri) and C(It%), respectively. 



2. Effective Hamiltonian results 

One of the main advantage of our approach is that 
we can study the stability of the variational ansatz 
wave function using the effective Hamiltonian method 
described in Sec. IIC. 17 Before showing our results, it 
should be noted first that one of the important quanti- 
ties which measure the quality of our approximation in 
the effective Hamiltonian approach is the effective con- 
stant K, K — 1 [determined as the result of Eq. (36)] 
meaning that there is no sign problem and the ground 
state of the effective Hamiltonian represents the exact 
ground state of H . As shown in the Tab. I, we found 
that the effective constant K is indeed very close to one 
for all the cases studied. This indicates that the num- 
ber of off-diagonal matrix elements, which cause the sign 
problem and are taken into account only approximately 
in H , represents only a very tiny fraction of the total 
number of matrix elements of H. Therefore, H is ex- 
pected to be rather close to H, and indeed it coincides 
with H in the large anisotropic limit J'/J — > 0, where 
there is no sign problem. 

Encouraged by this observation, we have calculated 
the spin-spin correlation functions C(l) and the dimer- 
dimer correlation functions D(l) using the FNE method, 
i.e., \ip} = |<$j ff ) in Eqs. (49) and (50), with the opti- 
mized |p-BCS) for \iPg) m H . Here IV'o ) is t ne exact 
ground state of H oS calculated numerically. The typi- 
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FIG. 8: System size dependence of (a) spin-spin correlation 
functions C(l) at the largest distance in the ri-direction (P a ) 
and (b) the dimer order parameter squared Pj (see the text) 
for clusters of L = / x /. The coupling studied is J'/J = 0.33. 
The variational Monte Carlo and effective Hamiltonian (FNE) 
results are denoted by squares and crosses, respectively. The 
straight dashed line in (a) is a guide to the eye. 
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FIG. 9: Dimer-dimer correlation functions D(f) with r = Ifi 
for J'/ J = 0.33 calculated using both variational Monte Carlo 
(squares) and effective Hamiltonian (crosses) methods. 



cal results for C{1) and D(T) are presented in Fig. 7 and 
Fig. 9, respectively. By comparing the variational re- 
sults with the FNE ones, it is evident that these corre- 
lation functions appear almost unchanged even at large 
distances, strongly suggesting that the projected BCS 
state |p-BCS) is very stable and accurate. This should 



be contrasted to the isotropic case (cf. Fig. 15), which 
will be discussed in the next subsection. 



A Systematic finite size scaling analysis of the order pa- 
rameters P s = C(l max Ti) and Pj calculated using FNE is 
also reported in Fig. 8 for J'/J = 0.33 and L up to 30x30. 
It is clearly seen in Fig. 8 (b) that Pj diminishes to zero 
in the limit L — > oo. Even though the FNE results for 
Ps shown in Fig. 8 (a) are almost the same as the ones 
for the variational calculations with |p-BCS), it is still 
difficult to rule out completely the possibility of a very 
small non-zero magnetic order parameter P s < 0.001 in 
the thermodynamic limit. Nevertheless, the fact that the 
spin-spin correlation functions calculated for L < 1000 
using the variational wave function |p-BCS) are almost 
identical to the ones calculated using the more accurate 
FNE approach (Fig. 7), strongly suggests that the mag- 
netic long range order is not likely to occur even for the 
more accurate FNE ground state |V'o ff )- 

From these results, we conclude that the ground state 
of the spin-1/2 antiferromagnetic Heisenberg model on 
the triangular lattice with J'/J < 0.6 - 0.7 (see also 
Sec. IV) is a spin liquid state with gapless spin excita- 
tions at the incommensurate momenta, described at least 
approximately by the present projected BCS wave func- 
tion |p-BCS). 

Finally, let us discuss the implication of our numerical 
study on the experiments for CS2CUCI4 where the esti- 
mated coupling is J'/J ~ 1/3. 11 Recent inelastic neu- 
tron scattering experiments for this material by Coldea 
et al 11 have revealed that the lowest spin excitation ap- 
pears at an incommensurate momentum (see Fig. 10), 
and the observed excitation spectrum consists of a broad 
incoherent continuum, at least, in the intermediate tem- 
perature region above the magnetic transition tempera- 
ture Tm- As we have shown in the case of the ID system 
(Sec. IIIB), the BCS elementary excitations naturally de- 
fine true spinous in our approach, and according to the 
gauge theory 48,54 they should behave as free fermions at 
low enough energy, namely close to the nodal points of 
the incommensurate momenta where = = (see 
Fig. 5). Therefore, the low energy spin excitations ob- 
served experimentally can be explained by a two-spinon 
broad continuum. At present, we cannot calculate di- 
rectly the dynamical spin correlation functions. How- 
ever, using the technique described in Sec. II D, we can 
calculate the lowest spin-1 excitation energy at each mo- 
mentum, and in Fig. 10 our results are compared with the 
experimental values estimated for CS2CUCI4 by Coldea et 
al. 11 As seen in Fig. 10, both results are in excellent agree- 
ment, considering that our calculated excitation spectra 
with L up to 30 x 30 appear rather well converged to 
the thermodynamic limit. This remarkable agreement 
strongly supports the presence of a spin liquid state in 
2D. 
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FIG. 10: Lowest triplet excitation energy as a function of 
momentum for J'/J = 0.33 and different cluster sizes L (in- 
dicated in the figures), calculated using the method intro- 
duced in Sec. II D. For comparison, the experimental val- 
ues (open circles) measured by inelastic neutron scattering 
on CS2CUCI4 11 are also plotted. Here Gi = An/^/i, and the 
experimentally estimated value of J = 0.37 meV is used. 11 
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FIG. 11: Nearest neighbor gap functions, A(ti) = A(t2) = 
A(r2 — fi) — A, which are consistent with the sign con- 
vention for the short range RVB state explained in App. D. 
The dashed (solid) bonds represent a positive (negative) sign 
[Eqs. (D8) and (D9)]. Notice that the unit cell is (2x1). The 
origin of the lattice (0,0) [Eq. (55)] is denoted by O. 



D. Isotropic triangular lattice with J' — J 

It is well known 3 that for the spin-1/2 antiferromagnet 
Heisenberg model on the spatially isotropic triangular 
lattice (J' — J), a faithful variational ansatz is the so- 
called short range RVB state |Vrvb) described by the 
following wave function: 



IVrvb) = 

{C} 



n im> 



(51) 



where the sum {C} runs over all possible nearest neigh- 
bor dimer covering of the triangular lattice (with equal 
weights and therefore implying that all the spatial sym- 
metries of the lattice are satisfied), whereas the product 
is over the corresponding nearest neighbor singlet pairs 
of spins located on each dimer (defining singlet valence 
bond configurations for a given dimer covering) between 
sites i and j, 



(52) 



sorted according to the lexicographic order (i < j). 

Indeed, very recent numerical calculations for the 
6x6 isotropic triangular antiferromagnet by the Lanczos 
method 21 have shown that the overlap between the short 
range RVB wave function |Vrvb) and the exact ground 
state |*o) is very large KV'rvbI^o)! 2 = 0.891, 55 and es- 
pecially the average sign, 



(S) 



\(x\* )\ 2 Sgn[(x\* )(x\^ 



RVB/ 



(53) 



is very close to the maximum value, namely, (S) — 0.971. 
These results clearly indicate that the phases of the ex- 
act ground state |^o) are very well described by the short 
range RVB ansatz state |Vrvb)- The values of the over- 
lap KV'rvbI^'o)! 2 and the average sign (S) are even much 
better than the ones corresponding to classical Neel or- 
dered wave functions considered previously, which also 
contain additional variational parameters. 4 It should be 
emphasized that an accurate description of the phases of 
the exact ground state |^o) by a simple variational wave 
function such as IV'RVb) is the essential ingredient for re- 
liable calculations based on the FN or FNE approach. 

Although IV'RVb) is a very good variational ansatz state 
for the ground state of the isotropic triangular antiferro- 
magnet, the present representation of this state [Eq. (51)] 
is rather difficult to handle and is not convenient to im- 
prove the state IVrvb) systematically by including, for 
instance, long-range valence bonds, because, within this 
representation, there is a "sign problem" even at the vari- 
ational level 3 . In order to treat more conveniently the 
short range RVB wave function |Vrvb) and its variants, 
we have derived in App. D an exact mapping between the 
short range RVB wave function |Vrvb) and the projected 
BCS state |p-BCS) on planar graphs, 56 namely, for most 
interesting lattice geometries including the triangular lat- 
tice, the square lattice, and the Kagome lattice. 

As shown in Apps. C and D, the short range RVB 
state |Vrvb) is equivalently described by the projected 
BCS wave function [p-BCS) with a special choice of the 
variational parameters: the only non-zero parameters are 
the chemical potential /1 and the nearest neighbor singlet 
gap functions Ajj with the appropriate relative phases 
shown in Fig. 11, and the limit of —fi 3> (Ajjl is as- 
sumed so that the gap function Ajj is proportional to 
the pairing function /y [Eqs. (C10)] considered in the 
App. D. 
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It is easily seen from Fig. 11 that the corresponding 
BCS Hamiltonian Hbcs is defined on a (2 x 1) unit 
cell, and thus Hbcs is not translation invariant. In 
fact, -ffscs is invariant under an elementary translation 
T2 : r — > f + T2 in the T2 direction, whereas it is not 
under an elementary translation T\ : f — > r + f\ in the 
T\ direction. Let us now show briefly that the transla- 
tion symmetry is recovered after the projection Vq, i.e., 
|p-BCS) is indeed translation invariant. To this end, one 
can combine the translation operation T~i with the SU(2) 
gauge transformation U : 



ct 

r,<7 



-ct 

r,cr 



(54) 



for r = m\T\ + U12T1 with mi odd. Under the compos- 
ite application of the transformations Ti^)U, the pro- 
jected BCS wave function |p-BCS) remains unchanged. 
Therefore, the |p-BCS) is translation invariant because 
the SU(2) gauge transformation U acts as an identity in 
the physical Hilbert space with singly occupied sites due 
to the projection Vg- 

Owing to this new, more convenient representation 
of the short range RVB state IV'rvb) by the projected 
BCS wave function |p-BCS), we can now calculate var- 
ious physical quantities using the standard variational 
Monte Carlo method. For example, the variational en- 
ergy £'(^rvb) of the isotropic triangular antiferromagnet 
is plotted in Fig. 12 for different system sizes. From 
the finite size scaling analysis of lattice sizes up to 
L = 30 x 30, we can easily provide a very accurate es- 
timate of the variational energy in the thermodynamic 
limit: E(^ RVB )/JL = -0.5123 ± 0.0001 for L -> 00 (see 
also Tab. II). 
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FIG. 12: Variational energy for the isotropic triangular lat- 
tice with J' I J = 1.0 calculated using the short-range RVB 
state I^rvb) (circles), the short-range RVB state with fi = 
(squares), the classical Neel state I^Ncci) (diamonds), 4 and 
the present projected BCS sate J^p-BCS) with spin Jastrow 
factor J s (triangles). 



Another important advantage of this projected BCS 
wave function representation is that it is easy to improve 



the variational ansatz state systematically. For instance, 
only by changing the chemical potential /i from a large 
value to zero, the variational energy is significantly im- 
proved as is shown in Fig. 12 and in Tab. II, where the 
variational state is denoted simply by IV'rvb) with /1 = 0. 
Notice that in this case the |p-BCS) is equivalent to 
a Gutzwiller projected free fermion state with nearest 
neighbor hoppings defined in a (2 x 1) unit cell. 57,58 

More generally, within the framework of the projected 
BCS wave function |p-BCS), we can easily extend the 
variational ansatz state to include long-range singlet va- 
lence bonds simply by adding non zero gap functions 
Aij's to the BCS Hamiltonian i?Bcs- 59 In order for 
|p-BCS) to preserve translation invariance, the pairing 
part of the BCS Hamiltonian is generalized in the follow- 
ing form: 



E 



H.c. 



(55) 



where the first sum runs over all lattice vectors f — 
r\T\ + r-iT-i (r\,ri\ integer), whereas the second sum 
Y^'t is for t m — TOifi + miTi (mi, m^. integer) with 
TO2 > or with mi > and m.2 = denoted by solid 
circles in Fig. 24 (a). It is readily shown that i? pa ir is 
invariant under T\(^IA and 7^. In order to minimize 
the number of variational parameters, we have assumed 
here that A(i m ) — A(lZ y t m ). Because of the phases 
of the gap functions in -ff pa i r (see Fig. 11 for the near- 
est neighbor gap functions), this BCS Hamiltonian is not 
guaranteed to be reflection invariant around the yz-plain 
(R y ). Nevertheless, as will be discussed later, our numer- 
ical calculations indicate empirically that the considered 
projected BCS state |p-BCS) becomes a fully symmetric 
spin liquid state only in the thermodynamic limit, a state 
being not only translation invariant but also reflection in- 
variant. It is interesting to notice that even though this 
I;? -BCS) is not fully symmetric on small lattice sizes such 
as a 6 x 6 cluster, a much better overlap | (p-BCSI^o) | 2 
as well as a much better average sign (S) are obtained for 
this |p-BCS) compared to the ones for a translation in- 
variant complex |p-BCS) ansatz. 19,20 In App. C, several 
peculiar and interesting properties of the corresponding 
unprojected BCS state |BCS) are derived analytically. 
For example, the spectrum of the BCS Hamiltonian 
has two branches [due to the (2x1) unit cell] and a small 
excitation gap, which vanishes when the long range gap 
functions are turned off (provided /i = 0). 

As mentioned above, it is extremely important to show 
that the projected BCS state |p-BCS), constructed from 
I BCS) with the paring interactions H pa i r , preserve all re- 
flection symmetries of the lattice in the thermodynamic 
limit, because otherwise some symmetry broken with fi- 
nite order parameter occurs, and this variational state is 
no longer a spin liquid. This symmetry property is eas- 



18 



Method (wave function) 


E/JL 




VMC (RVB) 


-0.5123 ± 0.0001 


0.0 


VMC (RVB with /x = 0) 


-0.5291 ± 0.0001 


0.0 


VMC (BCS+Neel) 60 


-0.532 ±0.001 


0.72 


VMC (present study) 


-0.5357 ±0.0001 


0.0 


FN 


-0.53989 ± 0.00003 


0.325 ± 0.006 


FNE 


-0.54187 ± 0.00006 


0.353 ± 0.007 


SW 


-0.538 ± 0.002* 


0.4774 


GFMCSR 4 


-0.545 ± 0.002* 


0.41 ±0.02 



TABLE II: The energy E and the magnetic moment m (di- 
vided by its maximum value mo = 1/2) estimated in the ther- 
modynamic limit for the spin- 1/2 antiferromagnetic Heisen- 
berg model on the isotropic triangular lattice (J — J'). In 
the first four rows are the variational Monte Carlo (VMC) es- 
timates for different wave functions: (from the top) the short- 
range RVB state IV'rvb), IV'Rvb) with fi — (see the text), 
the wave function studied recently by Weber, et af° in which 
a classical Neel order and a singlet pairing (with symmetry 
different from ours) are included, and the present wave func- 
tion J B \p-BCS). The fixed node (FN) and effective Hamilto- 
nian (FNE) results are then provided. For comparison, the 
corresponding values estimated by the linear spin wave the- 
ory (SW) and the Green function Monte Carlo method with 
stochastic reconfiguration (GFMCSR) 4 are also presented. 
Note that the energies computed by the last two methods 
(denoted by *) are not a rigorous upper bound of the exact 
ground state energy. 
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FIG. 13: System size dependence of the nearest neighbor spin- 
spin correlation functions C(f) [Eq. (49)] for different cluster 
sizes L = I x I up to 30 x 30 sites and for the three different 
directions of the triangular lattice (f — f\, fz, and T2 — f\). 
Here \ip) = \p-BCS) with A(n) = A(f 2 ) = A(f 2 - n) = 1 
[see in Fig. 11 and Eq. (55)], Uj = 0, and /i = —1. This state 
is a good but not optimal variational wave function for the 
spatially isotropic model (J = J'), and it is used only to show 
that wave functions of this type recover all spatial symmetries 
of the lattice in the thermodynamic limit. 



ily proved for the short range RVB case. Within open 
boundary conditions, the short range RVB state |V>rvb) 
[Eq. (51)] and the projected BCS state |p-BCS) with 



the gap functions defined as in Fig. 11, ti 



0, and 



fx — > — oo are exactly the same for any finite size clus- 
ters (see App. D). Therefore the symmetry of the state 
|p-BCS) is implied by the one of the the short range RVB 
state I^rvb)- For large clusters the boundary conditions 
cannot play a role, and hence the symmetry is approx- 
imately recovered even when periodic boundary condi- 
tions are used. We have checked numerically in Fig. 13 
that the reflection symmetry is very well preserved in 
the thermodynamic limit also for the general case where 
the pairing function fij is not restricted only to nearest 
neighbor bonds. However in this case we could not rigor- 
ously prove our empirical observation of this symmetry 
invariance because the equivalence of RVB wave func- 
tions and |p-BCS) no longer holds for long range pairing 
functions. 

In order to further improve our variational ansatz wave 
function, we also introduce a simple spin Jastrow factor 
Js into the projected BCS wave function described above, 
i.e., |p-BCS) -> Jsb-BCS), where 



Js = exp 



i,j=X 

(»<i) 



(56) 



and v(r)'s are variational parameters which are optimized 



using the SR method explained in Sec. II B (also see in 
App. A). Since the most important contributions to the 
variational energy are from the short range v(r) , s, in 
what follows we consider only terms with |r| < 3 and the 
v(r) , s for larger distances are set identically to zero. As 
will be discussed later, the inclusion of the spin Jastrow 
factor J s is also crucial for stable FN and FNE calcu- 
lations. Similarly, for the gap functions, only A(£ m )'s 
with |i m | < 3 are considered in the present variational 
wave function. After the SR optimization calculations, 
it is found that the optimized chemical potential /i is 
non zero and only nearest neighbor hopping terms with 
tij = 1 for all the directions are relevant for the isotropic 
case (J = J'). In Fig. 12 and in Tab. II, the calcu- 
lated variational energy is reported and compared with 
the results for other variational ansatz states. The im- 
provement of the present variational wave function com- 
pared with previously studied states is remarkable. Our 
variational energy is even sizably better than a very re- 
cent estimate reported in Ref. 60, in which a variational 
ansatz state studied includes both a classical Neel or- 
dered magnetic moment and a singlet pairing with sym- 
metry different from ours. To our knowledge, the present 
value for the energy per site in the thermodynamic limit, 
— 0.5357 J ± 0.0001J, is the lowest among the variational 
estimates reported so far. Of course, within the present 
ansatz, this value can be further improved by considering 
longer range terms in v(r) , s and A(£ m )'s. 
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It is worth mentioning here that in general the BCS 
excitation spectrum [Eq. (C6)] has a finite gap when 
the longer range gap functions are included in the BCS 
Hamiltonian as in the present case, and therefore, as op- 
posed to the spin liquid state discussed in Sec. Ill C, 
the spin liquid state described here by |p-BCS) generally 
shows a finite gap in spin excitations. It should be also 
noted that, although the variational ansatz state with 
the spin Jastrow factor J s breaks the SU(2) spin rota- 
tion symmetry, non-singular J s does not imply long range 
magnetic order. This is clearly seen in Fig. 14 where the 
spin structure factor 



-iq-r 



C(r) 



at q = (47r/3,0) is calculated for the optimized varia- 
tional wave function mentioned above. 
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FIG. 14: System size dependence of spin structure factor S(qj 
at q — Q* = (47r/3, 0) for the optimized variational state 
J7s|p~BCS) (see the text) for the spin-1/2 antiferromagnetic 
Heisenberg model on the isotropic triangular lattice of size 
L. Note that Q* corresponds to the momentum at which 
S(Q*)/L is finite for L — > oo when the state is classical Neel 
ordered with relative angle of 120° between the nearest neigh- 
bor spins on the different sublattices. 



Although the projected BCS wave function considered 
here is a spin liquid state with a finite spin gap and with- 
out any type of long range magnetic order (Fig. 14), 
when the FN or FNE method is applied with \tpa) = 
j7s|p~BCS), long range magnetic order appears with a 
finite magnetic moment in the z-direction (the effective 
Hamiltonian approach breaks the spin rotational invari- 
ance in this case and no detectable magnetic moment 
is observed in the other directions). Indeed, as shown in 
Fig. 15, the spin-spin correlation functions C(r) along the 



Ti direction are drastically increased for the FNE ground 
state \ipQ S ), compared with the ones for the spin liquid 
state J'slp-BCS). Furthermore, the oscillations observed 
in C(r) for the FNE ground state (Fig. 15) are consis- 
tent with the classical Neel order. These results strongly 
indicate that, in spite of its good variational energy, the 
spin liquid state is not stable against magnetic order for 
the isotropic case. 
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FIG. 15: Spin-spin correlation functions C(Zri) in the fi di- 
rection for the spin-1/2 antiferromagnetic Heisenberg model 
on the isotropic triangular lattice with L = 18 x 18 calcu- 
lated using both variational Monte Carlo (VMC) and effective 
Hamiltonian (FNE) methods. The variational wave function 
considered here is the projected BCS state ^|p-BCS) with 
spin Jastrow factor Js (see the text), and this wave function 
is used as the guiding function \^g) for the FNE calculation. 
The spin-spin correlation functions in the and r% — f\ di- 
rections are essentially the same as the ones presented here. 



The FN and FNE calculations of the spin-spin correla- 
tion functions were carried out by applying the " forward 
walking" technique developed in Ref. 33, which will be 
described briefly in the following. With this technique, 
the expectation value of any operator O diagonal in con- 
figuration space {x\, can be computed for the ground 
state of the effective Hamiltonian. This scheme removes 
completely the bias of the so-called mixed average esti- 



mate, 



61 



Omix = 



by a large "forward-walking" time r projection of the 
guiding function \iPg) m the above expression. More pre- 
cisely, the expectation value of O for the state |?/>Q ff ) is 
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computed through the equation 



cff\ 



WW! 



= lim 



{i> G \e- TH ^Q\i>f 
(^ G |e-^e ff |^ff) 



(57) 



where the simple relation lmv—^oo e~ THdl \ipo) oc \ipQ S ) is 
used. This scheme provides for a quantitative estimate 
of the magnetic moment (order parameter), as shown in 
Fig. 16, where one can clearly see that the values of the 
spin structure factor S(q) at the commensurate wave vec- 
tor q = Q* = (4/37r,0) considerably increase with clus- 
ter sizes and with the forward-walking projection time 
r, the r = value corresponding to the much less ac- 
curate mixed average estimate. In this figure, it is also 
apparent that a satisfactory convergence in t [Eq. (57)] 
is always obtained for the clusters studied. From the 
technical point of view, this calculation was made pos- 
sible due to the quality of our variational ansatz. For 
instance, without the inclusion of the spin Jastrow fac- 
tor J7s a much longer forward-walking time r is neces- 
sary to achieve a reasonable convergence, with an almost 
prohibitive computational effort to reduce the statistical 
errors to an acceptable value. 
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FIG. 16: Forward walking time evolution r (unit J -1 ) 
[Eq. (57)] of the spin structure factor S(q) at q = Q* = 
(47r/3, 0) for J'/J = 1.0. Here our best variational wave func- 
tion |p-BCS) with spin Jastrow factor J s (see the text) is used 
as the guiding function |V>g). For clarity, for each cluster, the 
FN results S(Q*)fn are referenced to the mixed-average es- 
timate .S^Q^mix, which is set to zero. 



Encouraged by these results, let us finally study the 
system size dependence of the spin structure factor S(q) 
at the commensurate wave vector q = Q*. The results 
calculated using both FN and FNE methods, S"(Q*)fn 
and 5'(Q*)fne, are presented in Fig .17. For each clus- 
ter, the FN and FNE Hamiltonians are constructed using 
the optimized variational state J s \p-BCS) as the guid- 
ing function IV'g)- I n this figure, to reduce the finite 
size effects for estimating the magnetic order parameter, 
the difference between the FN/FNE results S(Q* )fn /fne 



and the corresponding variational estimates 5(Q*)vmc is 
plotted. The difference should be extensive if there exists 
long range antiferromagnetic order, as the variational re- 
sults S(Q*)vmc/L for the spin liquid state J s \p-BCS) 
decreases to zero in the limit L — ► oo (Fig. 14). It is 
clearly observed in Fig .17 that the FN and FNE results 
for S(Q*)/L converge to finite values in the thermody- 
namic limit. 
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FIG. 17: System size dependence of the spin structure factor 
S(q) at q = Q* = (47r/3,0) for the spin-1/2 isotropic trian- 
gular antiferromagnet with J'/ J — 1.0 calculated using FN 
and FNE methods. Here our best variational wave function 
j7sjp - BCS) with spin Jastrow factor J s (see the text) is used 
as the guiding function \ipG)- In order to reduce finite size ef- 
fects, for each cluster, the variational result S(Q*)vmc/L for 
the spin liquid state ^Ip-BCS) (which is zero for L — > oo) is 
subtracted from the FN and FNE results S , (Q*)fn/fne/^- 



From these calculations we can estimate the magnetic 
moment quantitatively. To this purpose, it has to be 
taken into account that the long range magnetic order 
occurs in the z-direction because the FN Hamiltonian 
H FN and the FNE Hamiltonian H do not commute 
with the total spin operator and display magnetic or- 
der only in the z-direction. 62 Therefore, the magnetic 
structure factor is related to the magnetic order parame- 
ter (magnetic moment) m up to a simple factor, namely, 
5'(Q*)fn/fne/^ ~~ * m 2 /2. Our estimates of m are sum- 
marized in Tab. II along with the ones of the energy, Eq n 
and Eq S , in the thermodynamic limit. For comparison, 
the estimates calculated by other methods and for dif- 
ferent wave functions are also provided in Tab. II. It is 
remarkable that even though a spin liquid wave function 
is used here as a guiding wave function, a finite magnetic 
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moment m is obtained with the FN and FNE approaches. 
These results clearly indicate that the spin liquid is even- 
tually unstable in the isotropic model. It should be noted 
here that the magnetic moment estimated with approxi- 
mate techniques such as the FN and FNE methods is con- 
sidered as a lower bound because they are biased toward 
non magnetic order by the spin liquid guiding function. 

It is also important to notice in Fig. 17 and in Tab. II 
that the FNE approach, which is a better variational 
method than the standard FN method, estimates a siz- 
ably larger value of m compared with the FN result, much 
closer to previous estimates based on different guiding 
functions with explicit magnetic order. 4 In the previ- 
ous work, 4 though the energy was not rigorously vari- 
ational unlike in the present case, similar corrections to 
the guiding function were implemented. 63 It is therefore 
very interesting that, by using two different variational 
wave functions with (overestimated) or without (strongly 
underestimated) magnetic order, both with good varia- 
tional energy, and by applying very similar techniques, 
a finite order parameter m for the isotropic triangular 
lattice is obtained, rather independently of the guiding 
function used. This is a very important property of meth- 
ods, such as FN, FNE, and the previously introduced 
Green function Monte Carlo method with stochastic re- 
configuration, 4 which are all based on the numerical so- 
lution of an effective Hamiltonian. On the contrary, the 
simple variational approach does not lead to a reliable 
prediction for the magnetic moment m simply because 
very different variational ansatz states with or without 
a finite magnetic moment may have very similar energy 
(see e.g.. Fig. 12 and Tab. II) but opposite long distance 
behavior. This is an additional confirmation that the 
present approach, a systematic correction to the varia- 
tional ansatz, is very useful for understanding quantita- 
tive physical properties of strongly frustrated quantum 
systems, for which either numerically or analytically ex- 
act solutions arc not known. 



E. Nearly isotropic triangular lattice with J' < J: a 
spin liquid with a small spin gap 

In the previous subsection, we have shown a robust nu- 
merical evidence that a spin liquid state is not the ground 
state for the isotropic triangular antiferromagnet. How- 
ever, it is natural to expect that quantum fluctuations 
become enhanced by increasing the spatial anisotropy 
J/ J' > 1, and that the magnetically ordered state even- 
tually melts and a true spin liquid phase would emerge. 
In this subsection, we shall extend the spin liquid ansatz 
wave function discussed above away from the isotropic 
point, and study the stability of the spin liquid state 
with the FN or FNE method, as it was done successfully 
in the isotropic case. 



1. Stability against magnetic ordering 

As reported in previous studies, 5 a simple semiclassi- 
cal solution implies that the spin structure factor S(q) 
displays Bragg peaks at incommensurate momenta even 
slightly away from the isotropic case with J' / J ^ 1. It 
is confirmed by our variational approach, which provides 
a much better variational state compared to the classi- 
cal solution, that these peaks appear in S(q) and their 
locations in the Brillouin zone smoothly evolve from the 
commensurate momenta Q* = (An/ 3, 0) (and equivalent 
ones) to the incommensurate ones within our accessible 
finite size clusters. 

At the variational level with the same type of spin liq- 
uid wave functions considered for the isotropic case (see 
also Tab. Ill), the spin structure factor is finite for these 
incommensurate peaks, and therefore, as opposed to the 
classical solution, no long range magnetic order is im- 
plied. To study this property within the FN or FNE 
approach, we have to note that it is very difficult to per- 
form a finite size scaling analysis when incommensurate 
correlations are studied, a situation which is further com- 
plicated by the proximity to a possible phase transition 
from a magnetic to a non magnetic ground state, because 
the value of m, as we have shown, is rather small already 
in the isotropic case. In order to carry out a reliable fi- 
nite size scaling, in the following, we consider a sizable 
anisotropy within the hypothesis to be far enough from 
the critical point (which is unaccessible within the finite 
size clusters studied, L < 1000). Indeed, for J' / J = 0.7, 
we are able to successfully carry out the analysis of the 
FNE calculations, and the results are presented in Fig. 18. 
Here the guiding wave function |p-BCS) does not require 
the spin Jastrow factor J s to be accurate enough, and it 
consists of A(t m ) with \t m \ < 3 [Eq. (55)] as well as 
the hopping integrals and the chemical potential. All 
these parameters are optimized using the SR minimiza- 
tion method. As shown in Tab. Ill, it is found that the 
optimized chemical potential is non zero, and only the 
nearest neighbor hopping in the ri direction is relevant 
and can be set to unity (t^ =1). 

The excitation spectrum E^ of the corresponding BCS 
Hamiltonian are analytically calculated in App. C, from 
which several very interesting features are observed. If 
the gap functions A(t m ) are restricted to the near- 
est neighbors, the spectrum of the BCS Hamiltonian is 
generically gapless (provided fi — 0). Once the gap func- 
tions are non-zero for longer bonds, a tiny energy gain 
is obtained, and correspondingly a finite excitation gap 
in i?k is opened. However this gap is very small. In 
fact, it is found that the low energy BCS spectrum 
on finite size is almost indistinguishable from a gapless 
one, and that the lowest excitation gap in E^ is esti- 
mated as small as ~ 0.3% of 2W for J'/J = 0.7 where 
W is the maximum excitation energy of E^. 64 The rea- 
son for the appearance of a finite excitation gap in E^ 
is simply understood because the BCS Hamiltonian with 
broken translation symmetry has two sites per unit cell, 
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J'/J 


0.7 


0.8 




0.304(3) 


0.243(6) 


A(1,0) 


2.01(3) 


2.33(4) 


A(0,1) 


1.08(2) 


1.42(2) 


A(l,l) 


0.205(3) 


0.215(3) 


A(-l,2) 


0.002(3) 


0.001(3) 


A(2,0) 


-0.36(1) 


-0.32(2) 


A(0,2) 


0.01(1) 


-0.02(3) 


A(2,l) 


0.011(2) 


0.063(2) 


A(l,2) 


0.002(2) 


0.001(2) 


A(-2,3) 


0.001(1) 


0.0008(9) 




1.1886(3) 


1.2242(3) 


Evmc/JL 


-0.46467(3) 


-0.47840(3) 




-0.47051(2) 


-0.48521(2) 


-Efne/JX 


-0.47171(3) 


-0.48691(4) 



TABLE III: The optimized variational parameters of the 
projected BCS wave function for J'/J = 0.7 and 0.8 with 
L — 18 x 18. A(mi,m,2) is the singlet gap function A(i m ) 
for t m = miTi + m^n in Eq. (55), and a is the chemical 
potential. The remaining variational parameters are zero ex- 
cept for the nearest neighbor hopping in the fi direction, 
which is chosen to be one. The value of K, variational en- 
ergy Bvmc = E(p-BCS), FN (FNE) ground state energy 
£ FN = £o N (Efne = Ef) with \i> G ) = |p-BCS} are also 
presented. The number in parentheses is the statistical error 
bar corresponding to the last digit of the figure. 



and thus the BCS spectrum is in general gapped. It 
should be emphasized here that, as explained in the pre- 
vious subsection, after applying the projection Vg onto 
the GS | BCS) of this BCS Hamiltonian, the translation 
symmetry of the projected BCS state |p-BCS) is recov- 
ered. Therefore, a quite new, remarkable possibility to 
form a spin liquid state with a finite spin gap 65 and with 
a single spin per unit cell can be easily established here 
within the present variational framework without break- 
ing the translation symmetry. 

Let us now discuss the static magnetic properties for 
J'/J = 0.7. As seen in Fig. 18, the finite size scaling 
analysis of the FNE results for the spin structure factor 
clearly indicates a vanishing magnetic order parameter. 
We have also found that the FNE spin structure fac- 
tors as well as the FNE spin-spin correlation functions 
do not differ significantly from the ones computed for 
the variational spin liquid state [p-BCS) explained above 
(Tab. III). This should be contrasted to the isotropic case 
shown in Sec. HID, where significant differences were ob- 
served. These results strongly suggest that the spin liquid 
ansatz state described here by the projected BCS wave 
function |p-BCS) is stable within the present approach. 
Therefore, another spin liquid region different from the 
one discussed in Sec. Ill C might exist in the anisotropic 
triangular lattice with J'/J closer to one, the state of 
which can be described, at least approximately, by this 
proposed spin liquid wave function. 
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FIG. 18: System size dependence of the spin structure factor 
S(q) at q = Q* = (q*,0) for J'/J — 0.7 calculated using the 
FNE method. The incommensurate momenta q* as well as 
the corresponding value for the classical limit 5 are plotted in 
the upper panel. For clarity, the variational estimates of the 
structure factor S'(Q*)vmc are subtracted from the FNE re- 
sults 5 , (Q*)fne- The optimized projected BCS state |p-BCS) 
without spin Jastrow factor J s (see the text) is used for the 
variational calculations. The same wave function |p-BCS) is 
used as the guiding function \ipci) for the FNE calculations. 
The system sizes used are indicated in the figure. 



2. Stability against spontaneous dimerization 



Since the BCS Hamiltonian is defined with a (2 x 1) 
unit cell [Eq. (55)], it is natural to ask ourselves if a va- 
lence bond solid (not liquid) would be stabilized, which 
should also exhibit a finite spin excitation gap at the 
expense of a broken translation symmetry. In order to 
examine whether the valence bond solid is energetically 
more favored than the spin liquid state for the present 
anisotropic antifcrromagnetic Heisenberg model, it is suf- 
ficient to study the translation symmetry of the opti- 
mized projected state because the valence bond solid 
necessarily breaks this symmetry. This possibility can 
be easily considered within our approach by using a BCS 
Hamiltonian that is not invariant under the transforma- 
tion 7~i ^)U. In Eq. (55) it is assumed that the gap func- 
tion Aj j connecting sites i and j depends only on j — i up 
to the sign. To check the possible instability toward the 
valence bond solid, we have released these constraints 
for Aij but with remaining in the same (2x1) unit 
cell; for example, Aq ? x can be different from A^^, but 
Ao.fi = A2?!,3fi- We have then optimized independently 
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the first four Ay's at the shortest distances: 



' Ax 


— Ao^i, 


A 2 




' A 3 


= A ,f 2 , and 







(58) 



whereas all the other ones are kept fixed at the values ob- 
tained assuming the translation invariant ansatz (shown 
in Tab. III). Whenever the optimized parameters satisfy 
Ai = A2 and A3 = A4 within the statistical errors, the 
projected BCS state is a spin liquid. Otherwise the op- 
timized state is a valence bond solid simply because this 
state is no longer translation invariant. In order to opti- 
mize these variational parameters, we have used the SR 
minimization method described in Sec. II B. Our results 
are presented in Fig. 19, where the Monte Carlo evolution 
of these four parameters are plotted for J'/J = 0.7 and 
L = 18 x 18. It is clearly seen in Fig. 19 that, although 
the initial values of the parameters are far off form the 
symmetric condition (Ai = A2 and A3 = A4), after a 
few hundred SR iterations these parameters converge to 
the symmetric values (Tab. III). The inset of Fig. 19 
shows the Monte Carlo evolution of the corresponding 
energy as a function of the SR iterations (each iteration 
corresponds to a small variational Monte Carlo simula- 
tion with fixed variational parameters). After the first 
few SR iterations the energy appears to be trapped in a 
metastable state with a broken symmetry solution, but 
then after one hundred SR iterations the energy eventu- 
ally converges to a lower value corresponding to the spin 
liquid fully symmetric solution. These results, therefore, 
strongly indicate that the optimized state is translation 
invariant, and spontaneous dimerization is very unlikely 
to occur in this model because it is not stabilized even 
when the variational ansatz wave function allows to sta- 
bilize this kind of order. 

In principle, as shown in Ref. 52 for ID systems, also 
a translation invariant state can give rise to a sponta- 
neously dimerized order after applying the projection op- 
erator Vq on this state. In order to rule out this possible 
order, we have also calculated explicitly the dimer-dimer 
correlation functions D(r) defined in Eq. (50) for the vari- 
ational wave function described above (see also Tab. III). 
The results are compared with the ones computed using 
the FNE method. A typical example of the results is 
presented in Fig. 20 for J'/J = 0.7 and L = 30 x 30. 66 
As seen in Fig. 20, the two results do not show signif- 
icant differences, confirming that the spin liquid varia- 
tional ansatz appears very stable in this case. It should 
be emphasized once again that the situation is very differ- 
ent from the isotropic case where a magnetic instability 
was clearly detected with the FNE method (Sec. HID). 

The numerical calculations presented here strongly 
suggest that a new type of spin liquid discussed here is an 
appropriate ground state description of the spin-1/2 anti- 
ferromagnetic Heisenberg model on the anisotropic trian- 
gular lattice, at least in the region around J'/ J « 0.7-0.8 




FIG. 19: Monte Carlo evolution of the variational parameters 
defined in Eq. (58) as a function of SR iterations for the spin- 
1/2 antiferromagnetic Heisenberg model on the anisotropic 
triangular lattice with J'/J = 0.7 and L — 18 x 18. Here 
the SR minimization method explained in Sec. II B is used 
with St — 0.25/ J [Eq. (23)]. Inset: the corresponding energy 
evolution as a function of SR iterations. At each SR iteration, 
the energy is computed for the wave function with the fixed 
variational parameters given at the corresponding iteration in 
the main figure. 



(also see Sec. IV). In the present work, we have not at- 
tempted to determine the critical value of J' above which 
the incommensurate magnetically ordered state is stable, 
which was previously suggested in, e.g., Ref. 5. Indeed, 
there exists another more interesting phase boundary be- 
tween the two possible spin liquid states, the one pre- 
sented here and the one considered in Sec. IIIC, which 
will be discussed in the next section. 



IV. CONCLUSIONS AND FINAL REMARKS 

In this paper, using various quantum Monte Carlo 
techniques, we have studied the ground state phase di- 
agram as well as the low-lying spin excitations for the 
spin-1/2 antiferromagnetic Heisenberg model on the tri- 
angular lattice as a function of the spatially anisotropic 
coupling J'/J (Fig. 1). We have found numerical evi- 
dence for the presence of two different spin liquid states. 
The first spin liquid ("algebraic spin liquid") is stable 
for J'/J < 0.65 (see Fig. 21), and is characterized by 
gapless spin excitations, thus very similar to ID spin liq- 
uids. Conversely, the other spin liquid is rather a new 
type of spin liquid state, stable for 0.65 < J'/J < 0.8, 
and should show a small spin excitation gap. Starting 
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FIG. 20: Dimer-dimer correlation functions D(f) with r = Iti 
for the spin-1/2 triangular antiferromagnet with J'/J = 0.7 
and L = 30 x 30 calculated using both variational Monte 
Carlo (VMC) and FNE methods. The variational ansatz state 
considered here is the projected BCS state |p-BCS) (see text), 
and this state is used as the guiding function \ipa) for the FNE 
calculations. 



from the isotropic limit J' = J where the ground state is 
magnetically ordered (classical Neel ordered), quantum 
fluctuations increase strongly with decreasing the cou- 
pling J'/J down to zero. Therefore, the stability of a 
spin liquid in this region of the phase diagram is quite 
clear. 

The critical coupling J' c discriminating these two spin 
liquid phases can be determined in principle by com- 
paring the corresponding energy, and our best estimates 
of the energy as a function of J'/J are summarized in 
Fig. 21. From these results, it is concluded that the crit- 
ical coupling J'q/ J is about 0.65. Because of the limi- 
tation of currently available cluster sizes, our approach 
can not either determine precisely the critical coupling or 
describe accurately the nature of the transition. At the 
variational level, a first order transition occurs at a crit- 
ical point where the energy curves of the two spin liquid 
phases (as a function of J'/J), shown in Fig. 21, intersect 
with different slopes. In the same figure, it is also clear 
that the two slopes become very close within the FNE 
calculations, suggesting that the transition would eventu- 
ally turn to a conventional second order transition, when 
the quantum fluctuations are more accurately taken into 
account. 

Likewise, we have not tried to determine the critical 
coupling and to study the nature of the transition be- 
tween the spin liquid phase and the magnetically ordered 
phase with incommensurate magnetic order. The lat- 
ter phase should appear somewhere around J'/J > 0.8. 5 
However, it is very difficult to perform a finite size scal- 
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FIG. 21: Energy per site for the spin-1/2 antiferromagnetic 
Heisenberg model on the spatially anisotropic triangular lat- 
tice calculated using the variational Monte Carlo (VMC) and 
the effective Hamiltonian (FNE) methods. Here two types of 
spin liquid wave functions are considered, the one described 
in Sec. IIIC (denoted by "ID-type") and the other described 
in Sec. HID (denoted by "2D-type"). These wave functions 
are used as the guiding functions for the FNE calculations. 



ing analysis to determine the magnetically ordered in- 
commensurate phase simply because too large clusters 
are necessary for an accurate estimate of the magnetic 
moment. 

Our numerical results were obtained using the quan- 
tum variational Monte Carlo method as well as the lattice 
FN and FNE methods, which are essentially the Green 
function quantum Monte Carlo method with the fixed 
node approximation. The FNE ("effective Hamiltonian 
approach" ) is a further improved version of the standard 
FN method and developed here in Sec. II C. Although 
our results might be affected by this approximation in 
general, we have shown that the present methods provide 
very sensible results for the isotropic triangular antifer- 
romagnet with J' = J (Sec. HID), and the numerically 
exact results for the strongly anisotropic limit of the ID 
uncoupled chains with J' = (Sec. IIIB). Therefore, we 
have obtained rather reliable results for both limits of the 
phase diagram, and the same numerical tools have been 
applied also to the still controversial region of the phase 
diagram for < J'/J < 1. 

The quality of our approximation in the FN and FNE 
methods depends mostly on the choice of the guid- 
ing function |^>g)j f° r which we used the best varia- 
tional ansatz state, optimized using the SR minimization 
method (Sec. II B). The optimization of \ipG) is a cru - 
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cial ingredient of our approach, because the approximate 
ground state which we consider, i.e. the numerically ex- 
act ground state \ipQ S ) of the effective Hamiltonian H eS , 
can be computed with the restriction to have the same 
signs of IV'g)- Indeed, on small clusters for which numer- 
ically exact diagonalization of the systems can be done, 
the variational ansatz state described by an projected 
BCS wave function |p-BCS) provides a very good average 
sign (5*) [Eq. (53)] for both frustrated and non frustrated 
systems as discussed in Sec. HID and also in Refs. 4, 17, 
and 67. 

Within this approach, which was shown quite reli- 
able, we have found a surprisingly stable spin liquid 
("algebraic spin liquid") phase in the regime of large 
anisotropic couplings J/J' < 0.65 (Sec. IIIC). Our nu- 
merical calculations also indicate that this spin liquid 
state shows gapless, fractionalizcd spin excitations (Fig. 5 
and Fig. 6). 48,54 Therefore, we predict that this 2D alge- 
braic spin liquid state should show peculiar low energy 
properties similar to the ID systems. Although our con- 
clusion is based on numerical calculations for rather large 
clusters (up to 42 x 42 sites), to be fair, we cannot rule 
out a very weak instability toward symmetry-broken or- 
dered states, as predicted in the J' /J — > limit by us- 
ing the susceptibility criterion based on a random phase 
approximation (RPA). 68 Even in that study, the insta- 
bility occurs in an irrelevantly small temperature region, 
as also pointed out by the authors. 68 Moreover, it is not 
clear how reliable the RPA calculation is for finite values 
of J'/ J. 

We have also found another type of spin liquid phase in 
the region J'/ J close to the isotropic limit, i.e., for 0.65 < 
J'/ J < 0.8 (Sec. HIE). This rather new type of spin 
liquid state is characterized by a small spin excitation 
gap, and is described by the projected BCS state |p-BCS) 
with a gap function defined by Eq. (55) with a (2 x 1) 
unit cell. This spin liquid state is an extension of the 
conventional short range RVB state IV'rvb) [Eq. (51)], 
which has been considered before in the context of the 
isotropic triangular lattice 3 and is indeed a very good 
representation of the exact ground state of the isotropic 
triangular antiferromagnet for small clusters (Sec. HID 
and Rcf. 21). 

To extend the short range RVB state, we have con- 
structed an exact mapping between the short range RVB 
state and the projected BCS state |p-BCS) (Sec. HID 
and App. D). In doing so, the unit cell of the BCS 
Hamiltonian, defining the |BCS) state, is expanded to 
a (2 x 1) unit cell. This mapping is crucial in the present 
study because within this approach the short range RVB 
state can be easily extended and improved systematically 
by including hopping integrals, a finite chemical poten- 
tial, and most importantly long range gap functions in 
the BCS Hamiltonian, with no particular numerical ef- 
fort. An additional advantage of using the |p-BCS) rep- 
resentation is that it is easy to gain qualitative insight 
of the low-lying spin excitations from the corresponding 
BCS excitation spectrum E^- 48,54 For instance, E^ of 



the BCS Hamiltonian, corresponding to the short range 
RVB state, shows a finite excitation gap simply because 
— [i >• |Ajj| (Sec. HID and App. C), and therefore a 
finite spin gap is expected in the short range RVB state. 
This is indeed the correct property of the short range 
RVB state, because the presence of a finite spin gap and 
a very short correlation length have been established be- 
fore. 7 

It is worth mentioning a further remarkable property 
of this new type of spin liquid state described by this 
projected BCS wave function. Without the projection 
Vq, the BCS state |BCS) breaks translation and reflec- 
tion symmetries because the BCS Hamiltonian is defined 
with a (2 x 1) unit cell. As a consequence, the BCS 
excitation spectrum has a finite gap in general. How- 
ever, as discussed in Sec. HID, the translation symmetry 
is recovered after applying the projection operator onto 
|BCS). Therefore, within the projected BCS wave func- 
tions, we have discovered a peculiar way to open a finite 
spin gap without breaking the translation symmetry of 
the state. It is also important to emphasize that the re- 
flection symmetry of the projected state is also restored 
in the thermodynamic limit, as we have systematically 
tested numerically (Fig. 13), which however cannot be 
proved rigorously except for the limiting case of the sym- 
metric short range RVB state. 

Our finding of stable spin liquid phases in the spin- 1/2 
anisotropic triangular antiferromagnet is in good agree- 
ment with recent studies of the half-filled Hubbard model 
on the triangular lattice with spatial anisotropic hopping 
based on the Gutzwiller approximation 69 and the varia- 
tional Monte Carlo simulations. 70 However, it should be 
remarked here that the gap function in |p-BCS) which 
we found energetically favorable for the nearly isotropic 
case is rather different from the ones considered in the 
previous studies. 69 ' 70 Our spin liquid state is defined by 
a BCS Hamiltonian with a non translation invariant gap 
function, whereas the conventional ansatz state such as 
the ones considered in Refs. 69 and 70 is defined with 
a homogeneous gap function, which we found much less 
accurate close to the isotropic point (see Fig. 21). Al- 
though both the Hubbard model in the limit of large 
on-site repulsion U and the Heiscnberg model considered 
here should be the same, we just note that, with a vari- 
ational method sensitive only to the energy, it is much 
easier to find good variational wave functions for a low 
energy effective Heisenberg model rather than for the cor- 
responding Hubbard model, because the latter contains 
also the large energy scale U. 

In order to have better insight on the low-lying spin 
excitations of the spin liquid phases, we have directly 
calculated the spin one excitation dispersion with the 
method described in Sec. II D. This quantity is partic- 
ularly important since it can be compared directly with 
inelastic neutron scattering experiments. The detailed 
measurements on CS2CUCI4 by Coldea et al 11 are indeed 
available, for which a spin liquid like behavior has been 
observed. 11 It has been also reported that this material 
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can be described by the model Hamiltonian Eq. (1) with 
J'/J s» 1/3. 11 As shown in Fig. 10, our calculations were 
found to be in excellent agreement with the experiments 
with no fitting parameters. It is also clear from Fig. 10 
that the dispersion is 2D characteristic because for un- 
coupled chains (J' = 0) the dispersion should be sym- 
metric around the momenta (ir/2,0) and (w, 0). This ex- 
cellent agreement strongly supports the exitance of this 
type of 2D spin liquid state. The successful comparison 
also indicates that our numerical technique can compute 
accurately the excitation dispersion for the non-trivial 
strongly anisotropic antifcrromagnetic Heisenberg model 
on the triangular lattice. 

Encouraged by these results, we have also calculated 
the low-lying spin one dispersion for the other spin liquid 
phase which appears in the region of 0.65 < J' /J < 0.8. 
A typical results for J'/J = 0.8 is reported in Fig. 22. 
This coupling regime is relevant for the organic material 
k-(ET) 2 Cu2(CN) 3 for which a spin liquid like behavior 
has been also observed. 12 Since there is no experimental 
data available for the spin excitations of this material, 
our results shown in Fig. 22 provide a theoretical pre- 
diction for the spin excitation dispersion, which should 
be compared with the future neutron scattering experi- 
ments. It is interesting to notice that the dispersions for 
J'/J = 0.8 and for J'/J = 0.33 are very similar along 
the x-direction, whereas the difference becomes evident 
in other momentum regions, for example, along the y- 
direction. This feature should be also checked experi- 
mentally in the feature. 
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FIG. 22: Lowest triplet spin excitations as a function of mo- 
mentum for the spin- 1/2 antiferromagnetic Heisenberg model 
on the triangular lattice with J'/J = 0.8 and L — 18 x 18, 
calculated using the method introduced in Sec. II D. Here 
G 2 = 4tt/V~3 



In both spin liquid phases described by our variational 
ansatz, the spin excitation gap in the thermodynamic 
limit cannot be resolved by directly computing the excita- 
tion dispersions with the method introduced in Sec. II D 
as the present available system sizes are not large enough. 
According to the argument presented in the previous sec- 
tion, for J'/J ~ 0.8 a finite spin excitation gap is im- 
plied by the corresponding gap in the BCS excitation 



spectrum E-^ of the |p-BCS) ansatz. However, this gap 
should be very small because i) it is close to a transi- 
tion to an ordered state appearing for larger J'/J and 
ii) it is due to the rather small long range part of the 
gap functions Ajj's. In fact, we have found that the 
lowest excitation gap in E^ is as small as ~ 0.2% of 
2W for J'/J = 0.8 (where W is the maximum exci- 
tation energy in E^). 64 This small gap might also ex- 
plain an apparent finite spin susceptibility observed on 
k-(ET)2Cu2(CN)3, 12 which should eventually vanish ex- 
ponentially at very low temperatures with an activated 
behavior. 

It is also very interesting to note that the organic ma- 
terial k-(ET) 2 Cu2(CN) 3 , which shows a spin liquid like 
behavior under ambient pressure, 12 has been very re- 
cently found to become a superconductor under small ap- 
plied pressure (about 4-6 xl0 _1 GPa). 71,72 If we assume 
that the origin of the superconductivity is intimately re- 
lated to the spin liquid like nature of the phase observed 
next to the superconducting phase, the present study im- 
plies a rather unique, unexpected pairing symmetry of 
the superconductor with a non translation invariant pair- 
ing amplitude. This is a measurable effect because the 
projected BCS state |p-BCS) discussed in Sec. HIE is 
no longer translation invariant once the projection con- 
straint Vq is released or mobile carriers are introduced 
into the homogeneous spin liquid state. Therefore, we 
expect that this unconventional pairing formation can 
be probably observed by, e.g., scanning tunneling mi- 
croscopy experiments. 

Based on our numerical calculations and their compar- 
ison with experiments, we conclude that the spin liquid 
is a generic phase of quantum matter, and that, for its 
stability, it is not important to be very close to a phase 
transition, unlike the recently proposed scenario in which 
a spin liquid state appears only at a transition point be- 
tween a magnetically ordered phase and a spontaneously 
dimerized phase. 74 Within the projected BCS wave func- 
tion approach, spontaneous dimerization can be correctly 
described in quasi ID systems, but not in 2D unless the 
translation symmetry of the variational state is explicitly 
broken. 52 Even allowing this possible symmetry breaking, 
by adding appropriate symmetry breaking terms to the 
BCS Hamiltonian, we have not found a stable dimerized 
phase in the present 2D system. Instead, we have found 
rather stable spin liquid phases. Finally, our phase dia- 
gram of the spin-1/2 antifcrromagnetic Heisenberg model 
on the triangular lattice is summarized in Fig. 23, where 
possibly related materials are also indicated. In particu- 
lar, it would be very important to distinguish experimen- 
tally the two different spin liquid phases predicted here 
for different values of the coupling J'/J. 
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APPENDIX A: DETAILS OF THE NUMERICAL 
CALCULATIONS 

In this appendix, we will explain how to define ap- 
propriately and calculate efficiently the quantity Ok(x) 
[Eq. (12)] for each variational parameter ctk of the wave 
function |\&r aji .i ). Since the derivation is rather general 
and not restricted to the particular form of the wave func- 
tion used, we will here indicate the determinantal part 
of the wave function by |4>o), instead of |BCS), and thus 
l*{Qfc}) = ^g|*o)- Since it is generally found much more 
efficient, we will consider a wave function with definite 
number N of electrons, for which there are also impor- 
tant simplifications for deriving computationally conve- 
nient expressions for the quantities Ofc(x). 

To this purpose, we will use a particle-hole transfor- 
mation introduced by Yokoyama and Shiba, 75 which is 
defined by the following canonical transformation: 



4i = (-i) ii+i2 c i+L . 



(Al) 



Here c\ a (cj CT ) is an electron creation (annihilation) oper- 
ator at site ?i — i\T\ + i^Ti with spin <r(=|, J,), whereas 



C\ and Ci represent the new canonical operators corre- 
sponding to the original spin up (spin down) electrons 
for i < L (L < i < 2L). L is the total number of sites. 

After this particle-hole transformation (Al), it is clear 
that the BCS Hamiltonian [Eq. (3)] can be written as 



cs 



2L 

E3( 



^(hf; 



i,j=i 



i.j 



Cj 



(A2) 



where (fi r ( HF )) / are appropriate 2L x 2L matrix ele- 
ments, which can be straightforwardly computed from 
Eq. (3). Note that this matrix can be easily evaluated 
also for rather general Hartree-Fock, BCS Hamiltonians 
containing, e.g., non translation invariant terms, or sev- 
eral types of orders. The total number of particles A^ ph ) 
after the transformation (Al) is related to the total spin 
along the ^-direction S z — (iVj — Ny)/2 in the original 
representation through the following equation: 

2L L 

AT(ph) = £ C i Cl = J- ( c t TCiiT + Ci4 cl) =L + 2S Z . 

7=1 8=1 

Since S z is conserved in the original BCS Hamiltonian, 
the transformed Hamiltonian (A2) has a definite number 
jV(P h ) of particles, as anticipated. With this particle-hole 
transformation, it is possible to study the spin excitations 
of the BCS Hamiltonian with S z 7^ 0, as well as the 
ground state which belongs to the singlet S z — sector. 
Therefore, in the following, we will consider the general 
case with unrestricted A^ ph ). 

By using an appropriate unitary transformation, 
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a=l 
2L 



(A3) 



i=I 



Wbcs is diagonalized with a new set of quasiparticle op- 
erators {7! ,7a}: 



2L 



(A4) 



where, for convenience, the eigenvalues are sorted in as- 
cending order E\ < £2 < • • • < £2L- 

A natural choice of the Slater determinant part |<J?o) 
of the variational wave function is the ground state of 
Hbcs with A^( ph ) particles: 



(A5) 



where |0) is the vacuum (C/|0) = 0) of the Hilbert space 
after the particle-hole transformation. Thus the Slater 
determinant for a particle-hole transformed configuration 



\X)=C\C\ 



•CI 10) is 
(xl*n 



det [S] , 



(A6) 
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where S is a iV( ph ) x N^ ph ^ matrix, whose elements are 

(5) l>a = <0|^t|6) = (&) Jia> (A7) 

It's (1 < h < 2L) arc the "positions" of the par- 
ticles (I = 1, 2, ■ ■ ■ 7V( ph )) corresponding to the configura- 
tion |x), and a = {1,2,- • • ,M P ^}. 

Now let us consider small changes for etk — > a' fe = 
a.k + Jafc. Since 7Ybcs depends linearly on the 
perturbed system is described by 



"^bcs = "^bcs + E 5a k ■ V k 



(A8) 



k=l 



where V k = Eij=i {V (k) ) u C\Cj is a suitable opera- 
tor proportional to the chosen variational parameter otk ■ 
For instance, when the chemical potential is changed 
M — » M + ^ namely a term -^Eti E ff c [/'," 
is added to -Hbcs, the corresponding matrix element 
(VW) U is — &i,J for I < L and 5i,j for I > L, 
due to the particle-hole transformation (Al). Here 



is the Kronecker <5-function. 



When a small 
h.c. 



is con- 



is 



change in a pairing term 

sidered, the corresponding matrix clement (V^)jj 
(-l)h+h [Srfjj+L + 6i J+L Sj,i\. 

With the basis set formed by the quasiparticle opera- 
tors {tI, 7a}, the matrix is transformed as — > 
WV^U, and thus the first order correction to the state 
|*o) is easily computed as follows: 



l*o> 



where 
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\* )+O(6a 2 k ) 
(A9) 



From Eq. (All), it is now easy to evaluate the per- 
turbed Slater determinant (x| <!>[,), and thus an explicit 
expression for Ok{x) defined by Eqs. (12) and (13): 







: otherwise 



(A10) 

This expression is well defined as long as |3>o) is non 
degenerate, so that the denominator in the definition of 
is always non zero. This condition is rather generally 
satisfied for the BCS Hamiltonian. We can now readily 
express the state \&' ) in terms of the {C\C} basis set, 



l*o> = 



2L 



\<t>o)+0(5at) 



k=l 7,7=1 

(AH) 

where AfW = UQ^W. At each iteration of_thc SR 
minimization procedure described in Sec II B, has 
to be computed, because the operators {7 Q ,7 Q } change 
every time a new set of {ctk} is calculated. This can be 
done using four matrix-matrix multiplications. 
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Here Qji is a local single-particle "Green's function" , 

<*l<5)Al*°> 



Qji = 



E V) 



J a 



al ' 



(A13) 



which is computed and updated during the variational 
Monte Carlo iterations. Since the matrix does not 

depend on the configuration \x) and Qji is always known, 
only about L 2 operations are required to evaluate Ok (x) 
for each k and for each sampled configuration |x). 

Finally, we note that when the variational wave func- 
tion |*{ Qfc }) contains the spin Jastrow factor J s — 



exp 



t-ii,j=i 



(i<j) V *3 S * S j 



i.e., 



|* W} > = VgJs\*o), 
the calculation of Ok (x) corresponding to the variational 
parameter Vij is much simpler and straightforward, be- 
cause Ofc(x) is simply (x\SfSJ\x) in this case. 



APPENDIX B: MARSHALL SIGN RULE AND 
PROJECTED BCS WAVE FUNCTIONS 

In this appendix, we will show that a BCS wave func- 
tion projected onto the subspace of singly occupied sites 
satisfies the Marshall sign rule provided the correspond- 
ing BCS Hamiltonian 



: for r\ > iV(P h ) and v < 



Hbcs = E 

3,1 



(E c ^ c ^ + ( a j.' c 5t' 



Let, + h.c. 



satisfies the following conditions: 



(Bl) 



j — : for j and I on the same sublattice. 

(B2) 

Here c\ a [ci a ) is an electron creation (annihilation) oper- 
ator at site Yi — {i Xl i y ) with spin g{=] 1 [), and tj.k and 
Aj^k are assumed symmetric under j <-> k interchange 
and real. It is also assumed that the ground state of the 
BCS Hamiltonian |BCS) is unique, namely with a finite 
size gap to the first excitation, a condition which can be 
generally met for non trivial values of Ajj and Uj. 

In what follows, we will show that the projected BCS 
wave function |p-BCS) = "Pel BCS) constructed from the 
BCS state above satisfies the Marshall sign rule, namely, 



Sgn[(x|BCS)] = (-l) 



N A 



(B3) 
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where is the number of down spins on one of the two 
sublattices, and \x) is an electron configuration with no 
doubly occupied sites, which can be written as 



<zi=<oin« 



i=l 

with the Ci ai factors being ordered from left to right ac- 
cording to the increasing index i (at is the spin of the 
electron on site i). Here it is assumed that the number 
of sites (L) is even. 

First, it is convenient to perform the following particle- 
hole transformation for the down spins only: 



(B4) 



With this transformation, the BCS Hamiltonian is trans- 
formed to a standard one-body Hamiltonian commuting 
with the total number of particles, i.e., 



Hbcs = ^ 



3,1 L 



^,1 



E4a 



]d\di al 2iS i 



v 

(7(7' 



,(B5) 



where "z" denotes the imaginary unit, S y is the y- 
component of the spin matrix, and condition (B2) is used. 
On the other hand, after this transformation (B4), the 
basis of no doubly occupied sites turns into 



I) 



|o>, 



namely, the sites are either doubly occupied or empty 
in the new representation, and all configurations {£} are 
now defined by the L/2 positions of the doubly occupied 
sites: 

(x\ = (0|dR 1 j.dR 1 tdR 2 j.dR 2 i- ■■ ■dR L/2 |dR Il/2 T) 

where R; € {ri,r2,--- , r^}. Because of the addi- 
tional phase in Eq. (B4), (x\ is related to (x\ by (x\ = 
(— 1) Na (x\. Therefore, using the definition of the Mar- 
shall sign (B3), a state with the Marshall sign rule in the 
original basis {x} is equivalent to a bosonic state in the 
particle- hole transformed basis {x}: for each configura- 
tion \x), the wave function should be always positive (or 
always negative). 

In order to show this bosonic rule for the projected 
BCS state considered, it should be noticed that the 
new basis {x} is invariant under a global rotation of 
the spins which transforms the spin matrix S v , into 



<r<W/2: 
d 3\ 



-( 1 

V2 V -* 



a 3l 



(B6) 



Indeed, dj^dji is transformed into aj^aj^. It is also eas- 
ily shown that this transformation (B6) factorizes the 
Hamiltonian (B5) as 



Hbcs = ^2 



E 



mi),," 



(B7) 



where = Hf, and being appropriate one- 
body Hamiltonian matrices whose eigenstates are </> a (r) 
and 0„(r), respectively (a = 1, 2 • • • , L). Therefore, the 
ground state |BCS) of Hbgs once computed in the basis 
of doubly occupied and empty sites {x} reads 



<x|BCS) = |Dct (S)\ 2 > 0, 



(B8) 



where S is a (L x L) matrix with (S 1 ), q = <^ a (R;). This 
proves the statement given at the beginning of this ap- 
pendix. 

Finally, we note that the condition (B2) is satisfied 
whenever the BCS Hamiltonian (Bl) is invariant under 
a particle-hole transformation: 



(-1) 



v r ■ 



(B9) 



provided tj t k — tkj and Aj,k is real. Eq. (43) in Sec. Ill A 
follows readily from this particle-hole symmetric condi- 
tion. 



APPENDIX C: BCS WAVE FUNCTION WITH A 
BROKEN TRANSLATION SYMMETRY 

In this appendix, several important properties are dis- 
cussed on the ground state |BCS) of the BCS Hamil- 
tonian which breaks primitive lattice translation sym- 
metry by extending the unit cell to (2 x 1). The re- 
sulting state |BCS) is used to build a projected BCS 
states |p-BCS) = "Pg|BCS) for the isotropic and nearly 
isotropic triangular systems discussed in Sec. Ill D and 
Sec. HIE, respectively. 

Let us start with defining the BCS Hamiltonian Hbcs 
introduced in Sec. II A in a slightly different way. A part 
of Hbcs which is invariant under any primitive lattice 
translations on the triangular lattice (Fig. 1) can be gen- 
erally described by the following Hamiltonian: 



MO) 

BCS 



E 



(Cl) 



where the first sum Y]-, runs over all lattice vectors f — 
T\T\-\-r2T2 {r\,r2~- integer), whereas the second sum X^F m 
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is for t m — m\T\ + m 2 r 2 (mi, m 2 : integer) with m 2 > 
or with mi > and m 2 — as denoted by solid circles in 
Fig. 24 (a). Now let us add to an additional pairing 





-o-o-o-^-o-o-o 
-o-o-o-o-o-o-o 

-o-o-o-i-o-o-o 



FIG. 24: (a): Region of (mi, 7712) (denoted by solid circles) 
considered for the sum over t m = m\T\ + 7712 T2. Each cir- 
cle corresponds to a pair of integers (7711,7712). (b): The first 
Brillouin zone (hexagon drown by dashed lines) for the trian- 
gular lattice shown in Fig. 1. The corresponding reciprocal 
lattice vectors are gi = 27r(I, — -^=) and 52 = 27r(0, -^j)- The 
reduced Brillouin zone (tilted rectangle) is also given by solid 
lines. Several symmetric points are F (0,0), P: (|7r, 0), Q: 



P':(i*,» 



and Q': (0, ^ tt) 



term Vp a i r which breaks the underlying lattice translation 
symmetry: 



E 



(C2) 



+H.c. 



Note that the unit cell of the total BCS Hamiltonian 
Vpair is now extended to (2 x 1). In the 
Aj- , Aj* , and \x are all assumed to be 
real. 



H 



BCS 



BCS 



following, tj- 



Let us first explore the excitation spectrum of Hbcs- 
For this purpose, it is convenient to Fourier transform the 
BCS Hamiltonian to the momentum space with the recip- 
rocal lattice vectors, gi — 27r(l, — ^) and g 2 — 2-7r(0, -^) 
[see Fig. 24 (b)]. After the Fourier transformation with 
^ra = -JzY<k e ~ lk ' Fc \ a ' Hbcs can be conveniently de- 
scribed by the following (4 x 4) matrix form: 



#bcs - ( 4f 4 



'fe+QT 



c -fe|> c -(k+Q)i 



( a(fc) Au(fc) A 12 (fc)\ / 

£ 2 (fc) A 12 (fc)* A 22 (fc) 

An(fe) A 12 (fc) -Ci(fc) 

VA 12 (fc)* A 22 (fc) -Uk)J \ c '-(Uq)iJ 



C fe+QT 



(C3) 



where 



£i(fc) = -2^ t t - m cos (k ■ t„^j - 11 
&{k) = -2^2 i t - m cos \(k + O) ■ i 



and 



' An(fc) = 2^' A t* m cos ■ t™) 



A 22 (fc) = 2^A t - m cos 



+ Q ) • t 



(C4) 



(C5) 



Ai 2 (fc) = ^A t - m (e<^)' f ™ + e 
with Q = gi/2 so that Q ■ r — nri, The primed sum 



r 



runs over the reduced Brillouin zone shown in Fig. 24 
(b), or equivalently, e.g., for parallelogram lattice of 
L = Li x L 2 (£1, L 2 : even) with k = ^g ± + &g 2 , 
the primed sum is taken over ki = 0, 1, • • • , L±/2 — 1 and 

fc 2 = 0,l,--- ,£ 2 -i. 

The excitation spectrum of the BCS Bogoliubov mode 
for this BCS Hamiltonian is now easily calculated by 
diagonalizing the (4 x 4) Hamiltonian matrix given in 
Eq. (C3). It turned out that the excitation spectrum E^ 

is doubly degenerate at each momentum k and is given 



by Er = E^, where 



Q a (k)±Q 1 {k) 



1/2 



(C6) 



and 
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Qo(k) = 
Qi(fc) = 



ti(k) 2 + Ukf + An(fc) 2 + A 22 (fc) 2 + 2 A 12 (fc) 

_ \2 



(Hi(k) 2 - + An(fc) 2 - A 22 (fc) 2 )' + 4 |A 12 (fc)N (a(fc) - 6(fc)) 2 + (An(fc) + A 22 (fc) 



- \ 2 



1/2 



r 



From this expression, it is clear that generally E^ has 
a finite excitation gap except for some special cases. 
So far no assumption has been made for the values of 



real gap functions Al- 



and 



As will be discussed 



in App. D, in order to be consistent with the phase rule 
[Eqs. (D8) and (D9)] and for Hbcs to be invariant under 
the transformation T^W (Sec. HID), the gap functions 
are subject to the following conditions: 



A t - m and A^ 



0, for m 2 odd 
0, for m 2 even. 



(C7) 



With this condition, A(t m ) given in Eq. (55) is equivalent 



to A t - and Af defined here. 

Let us now consider one special but important case 
for which all hopping terms are zero and An (ft) = 

— A 22 (fc). In this case, Q\{k) = 0, and therefore the 
excitation spectrum for the Bogoliubov mode is simply 



fJ 2 + A 11 (k) 2 + A 12 (fc) 



1 1/2 



(C8) 



After tedious but straight forward calculations, the 
ground state |BCS) of the BCS Hamiltonian Hqcs f° r 
this special case can be calculated analytically: 



| BCS) = exp 



E 

k 



k] C -(k+Q)l 



'hi -(k+Q)tJ) 



|0) (C9) 



r 



where 



fn(k) = 

/l2(fc) = 



-/22(fc) = 

Aia(fc) 



An (A) 



Ek 



(CIO) 



and A i2 (fc) = A 12 (—k) is used. This state should 
be compared with the more conventional one given by 
Eq. (9) where there is only one site per unit cell. Fi- 
nally, let us consider two more specific cases separately. 



Case (1): fi = and only nearest neighbor gap func- 
tions arc finite, i.e., A^ = A? 2 = A^ 2 _^ = A (see 
Fig. 11). Because An (fc) = -A 22 (fc) = 2Acos(fc-ri) and 
Ai 2 (fc) = 2A cos(fc • f 2 ) — i sin (k ■ (f 2 - fifj , has 

gapless excitations at fc = (/i/4±<? 2 /4 and — gi/4± g 2 /4. 
Case (2): /j, — > — oo. Since /n(fc) cx — An(fc) and 
/i 2 (fc) cx — Ai 2 (fc) in this limit, i.e., a pairing function 
is proportional to a gap function, a projected BCS wave 
function |p-BCS) built from this BCS state becomes 



J 



lim |p-BCS) = V G 

fl^ — OC 



E E' (Ac + (-ir Arj ( c t T c^ i - c t iC t +fmT ) 



L/2 



|0). 



(Cll) 



r 



From this result, it it clear that the projected BCS state 
discussed here indeed includes a short range RVB wave 
function IV'rvb) since if A^ = Ay 2 = A^ 2 _^ (all other 



gap functions are zero) is chosen, the above wave function 
(Cll) consists of the pairing functions with exactly the 
same phase as in Eqs. (D8) and (D9). More details of 
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this relation will be found in the next appendix. 



APPENDIX D: PFAFFIAN AND PROJECTED 
BCS WAVE FUNCTIONS 

We start from the definition of a Pfaffian in terms of 
an antisymmetric 2N x 2N matrix / where fij = —fj,i- 
The Pfaffian of matrix / is 

N 

p[f}= E ( D1 ) 

(»l<jl)i(*2<j2),— ,(«JV<JJV) fe=l 

and ii<«2'"<*JV 

where the sum runs over all possible covering of indices 
{{h,ji), {h, h), ■■■ , (inJn)} such that i k < j k and i x < 
12 < h < ■ ■ ■ < in, being p the parity of the permutation 



of the 2N indices: 

f 1 2 3 4 • • • 2N — 1 2N \ 
V h ji «2 h ■ ■■ *n 3n ) 

The most important relation known for the Pfaffian is 

that (P[f]) = Dct[/], which however will not be used 
in the following. 

Now let us suppose that the indices 1, 2, 3, 2N la- 
bel the positions of a lattice (not necessarily one dimen- 
sional). Then each covering of the indices in the Pfaffian 
is interpreted as a particular dimer configuration in which 
for example a spin singlet pair located at sites (ik,jk) is 
assigned. 

We can now define two spin wave functions in terms 
of these dimer coverings. The first wave function is ex- 
panded in the well known valence bond basis: 76 



|RVB) 



E 



i-ir 



(»l<jl),( i 2<J2),--- ,(lN<jN) 

and i\<i^---<iN 



N 



Lfc=l 



(D2) 



where S i is the spin- 1/2 lowering operator at site i, \F) 
is the ferromagnetic state defined by 



"It ' 1 



.t 

"2,-p 



-2N 



T |o), 



(D3) 



and |0) is the electron vacuum. Notice that for conve- 
nience we have included in the definition of the wave 
function the same permutation sign (— l) p appearing in 
the Pfaffian (Dl). The second wave function is the pro- 
jected BCS state: 



|p-BCS) =-p G exp 



27V 

E 

(i<3) 



hi ( c 1t c ],4 c L c ],t) 



|0> 



(D4) 



where Vq is the Gutzwiller projection operator onto 
singly occupied sites. In both wave functions, we have 
used the upper triangular part of the antisymmetric ma- 
trix /, whereas the tedious permutation sign (— l) p is 
present only in the wave function |RVB). 

We shall now show that the two wave functions 
|p-BCS) and |RVB) are actually the same for any f and 
on any lattice with all possible boundary conditions. 



Proof. Because the projection Vg forbids double occu- 
pancy and thus singlet bonds sharing the common site, 
the projected BCS wave function |p-BCS) can be ex- 
panded for all dimer coverings as follows: 



IP-BCS) = 



2N 

E 

. (*< j) 



'3,1 



c 5.t) 



JV 



10) 



E 

L<Jl) ) ( l 2<j2) ) - 

and i\ <«2' 

i 



>(»N<jiv) fc= 



N 



ikA°jk,l 



tk,i 



, T ) |o>, 



_[D5) 



where in the last formula the constant 1/N\ cancels out the ferromagnetic state in the wave function |RVB) and 
the factor due to the ordering of the indices i\ < i<i, • • • < making the necessary p* fermion permutations for each 
ifq. On the other hand, by substituting the definition of dimer covering, we arrive at the following expression: 
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|RVB) 



(»l<jl)i(»2<j2),-" ,(iN<jw) 

and ii<l2---<iN 



(-l)f(-lf 



JV 



|0). 



(D6) 



By simple inspection, it is readily realized that the 
fcrmion sign (— l) p is exactly the one (— l) p defining 
the permutation sign of the covering of indices {(ik,jk)}, 
and that Atl°> = c^cj^ -cj^cJjO). 

Thus it is proved that two wave functions |RVB) and 
|p-BCS) are exactly the same. 

In the remaining of this appendix, we will show that 
even for the triangular lattice geometry, the so-called 
short range RVB wave function IV'rvb) [Eq. (51)] can 
be described by the projected BCS wave function (D4) 
with a particular choice of the matrix /. 



for k — (0,0), whereas 



fk, : 



1 for j = (2,0) 
-1 fori = (1,1) 
-1 for j = (0,1) 



(D9) 



for k = (1, 0). All the other values of fkj are obtained by 
translation of 2f\ = (2,0) and/or r 2 = (1,0). Therefore, 
the unit cell of fk j is (2 x 1) on the triangular lattice (see 
Fig. 11). 



Consequences of the theorem: short range RVB 
wave function 



The short range RVB wave function I^rvb) [Eq. (51)] 
is defined with the same weight for all nearest neighbor 
valence bonds. Therefore, if we can find the matrix / in 
Eq. (D2) in such a way that the sign of the permutation 
for each dimer covering is exactly canceled, i.e., 



N 



(-i) p n^* =i > 



(D7) 



k=i 



then a relation is established between RVB wave func- 
tions and projected BCS wave functions. The condition 
(D7) is highly non trivial and difficult to satisfy as the 
number of dimcr coverings is exponentially large and the 
entries of the matrix are few in comparison. Neverthe- 
less, this problem for the case of nearest neighbor dimcr 
covering was solved in all planar graphs. 77 

By applying these old results, 77 we can generalize the 
Read-Chakraborty relation between the short range RVB 
and projected BCS wave functions on the square lat- 
tice 78 for the triangular lattice case. Namely, by using 
the known matrix / reported in Rcf. 77, it is possible to 
satisfy the condition (D7) even for the triangular case. 7 
Here we use open boundary conditions for a lattice oflxl 
sites (I being multiple of 6) where site i is ordered lexi- 
cographically: i = lm 2 + mi = (mi, 7712). The Cartesian 
coordinates of the triangular lattice (Fig. 1) are thus 

n = m\T\ + m 2 f 2 = (mi + m 2 /2, V3 m 2 /2). 

The pairing functions fk,j which meet the condition (D7) 
reads: 



1 for j = (1,0) 
fkj = \ 1 for j = (0, 1) 
1 for 3 = (-1,1) 



(D8) 



Complex representation 



We can multiply by the imaginary unit i all the a 
for the odd f\ components of j since the resulting wave 
function is equivalent to the original one [Eq. (D4)] in 
the presence of the projection operator Vg (apart from 
an overall phase). The obtained new complex matrix / 
consists of: 



fk,. 



for k = (0,0), whereas 




(1,0) 

(0,1) 

(—1, 1), the diagonal bond 



fk, : 



i fori = (2,0) 
1 for i = (1,1) 

— i for j = (0, 1), the other diagonal bond 



for k — (1,0). All the other values of fkj are obtained 
by translation of 2ti = (2,0) and/or f 2 = (1,0). 

Notice that if the diagonal bonds are eliminated, 
we end up with again a (1 x 1) unit cell with /k oc 
cos(fc x ) — icos(ky) in Fourier space, thus recovering the 
Read-Chakraborthy result for the square lattice with 
f t = (mi, m 2 ). 78 

It is important to emphasize that for the triangular lat- 
tice geometry a (2 x 1) unit cell is inevitable even with the 
complex representation (as the diagonal bonds acquire 
different signs). Thus, we conclude that the translation 
symmetry in fkj has to be broken for a good projected 
BCS wave function in the triangular case. It is interest- 
ing that two sites per a unit cell is enough, in contrast to 
the classical Neel order state with 120° degrees of mutual 
spin orientations which contains instead three sites per 
unit cell. 
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3. Periodic boundary conditions 

For the short range RVB state with periodic bound- 
ary conditions, it is known that condition (D7) can not 
be satisfied by a single |p-BCS), but by four, all ob- 
tained with the possible choices of periodic and antiperi- 



odic boundary conditions in t\ and t-i directions. 77 The 
proof can be applied for each of such |p-BCS), showing 
that these four projected BCS wave functions have to be 
used to exactly match the short range RVB wave function 
with periodic boundary conditions. 
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